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ABSTRACT 


Atomic  Clocks  provide  "stable"  signals  that  are  mainly  used  to  generate 
time  scales  and  to  measure  differences  of  time  between  events.  Each  atomic  clock  can 
individually  be  characterized  according  to  the  stability  of  the  scale  it  produces. 

Due  to  the  stochastic  behavior  of  each  clock,  usually  clock  ensembles  are  used  to 
build  more  stable  time  scales.  This  process  requires  basically  two  steps.  First  it  is 
required  to  individually  characterize  each  time  source  to  identify  the  particular  noise 
components  present.  Second,  a  measure  of  perfonnance  is  required  in  order  to  derive  an 
algorithm  based  on  it  to  properly  "weigh"  this  particular  clock  in  the  process  of  creating  a 
combined  scale. 

In  this  thesis  both  problems  are  faced  using  an  operations  research  approach:  each 
clock  is  modeled,  analyzed  and  characterized  by  a  time-dependent  measure  of 
performance  related  to  its  intrinsic  stability,  and  optimally  combined  to  produce  a  more 
stable  combined  time  scale.  The  optimality  criterion  is  directly  related  to  the  spectral 
characteristics  of  the  noise  sources  present. 

Real  laboratory  data  provided  by  the  “Real  Observatorio  de  la  Armada”,  ROA 
(Royal  Observatory  of  the  Spanish  Navy)  is  used  to  tune  the  final  algorithm  to  the  time 
standards  present  in  the  ROA. 
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EXECUTIVE  SUMMARY 


Atomic  Clocks  provide  "stable"  signals  that  are  mainly  used  to  generate 
time  scales  and  to  measure  differences  of  time  between  events.  Each  atomic  clock  can 
individually  be  characterized  according  to  the  stability  of  the  scale  it  produces. 


Due  to  the  stochastic  behavior  of  each  clock,  usually  clock  ensembles  are  used  to 
build  more  stable  time  scales.  This  process  requires  basically  two  steps.  First  it  is 
required  to  individually  characterize  each  time  source  to  identify  the  particular  noise 
components  present.  Second,  a  measure  of  perfonnance  is  required  in  order  to  derive  an 
algorithm  based  on  it  to  properly  "weigh"  this  particular  clock  in  the  process  of  creating  a 
combined  scale. 


In  this  thesis  both  problems  are  faced  using  an  operations  research  approach:  each 
clock  is  modeled,  analyzed  and  characterized  by  a  time-dependent  measure  of 
performance  related  to  its  intrinsic  stability,  and  optimally  combined  to  produce  a  more 
stable  combined  time  scale.  The  optimality  criterion  is  directly  related  to  the  spectral 
characteristics  of  the  noise  sources  present. 


The  thesis  has  been  divided  in  four  main  sections  with  the  following  contents: 


Section  one  is  a  general  introduction  of  the  problem  of  characterization  of  clocks 
and  precision  oscillators.  The  theory  of  analysis  of  stability  and  the  physical  magnitudes 
involved  in  this  process  are  also  presented. 


xv 


The  different  kinds  of  noise  processes  that  may  be  present  in  precision  oscillators 
are  outlined  and  their  intrinsic  characteristics  in  the  frequency  and  time  domains  are  also 
introduced  as  background  theory. 


Special  attention  is  given  to  the  main  statistic  used  in  frequency  stability:  the 
Allan  variance.  The  need  of  a  new  estimator  to  analyze  the  time  series  of  measurement 
data  is  shown  and  illustrated. 


In  section  two,  the  complete  characterization  problem  is  presented  as  a  nonlinear 
optimization  problem  subjected  to  nonlinear  constraints.  The  computation  of  weights  to 
assign  to  each  clock  with  the  objective  of  finding  a  composite  time  scale  is  also 
addressed. 


Classical  approaches  are  presented  and  a  new  approach  is  outlined. 


In  chapter  three,  simulation  is  used  to  test  and  compare  the  various  approaches. 

Four  different  tests  are  carried  out: 

•  A  single  point  simulation  for  a  hypothetical  case  that  has  been  analyzed 
previously 

•  Random  generation  of  covariance  matrices 

•  A  test  of  sensitivity  carried  out  on  the  new  proposal  to  check  if  small 
changes  in  the  input  parameters  lead  to  coherent  results 

•  Finally,  a  complete  simulation  involving  the  generation  of  time  difference 
data  under  the  presence  of  noise  is  perfonned.  The  importance  of  the 
length  of  the  time  series  is  addressed. 


xvi 


Finally,  in  section  four,  the  method  is  applied  to  real  data  provided  by  the  Royal 
Observatory  of  the  Spanish  Navy.  A  complete  characterization  of  an  ensemble  of  five 
atomic  clocks  is  performed. 


XVII 
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I.  INTRODUCTION.  CHARACTERIZATION  OF  CLOCKS  AND 

OSCILLATORS 


A.  INTRODUCTION 


Precision  oscillators  play  an  important  role  in  many  aspects  of  our  daily  life. 
High-speed  communications,  space  tracking,  navigation  and  many  other  important 
applications  are  based  on  a  precise  and  accurate  synchronization  among  the  different 
devices  and  the  equipment  that  is  involved  in  such  activities.  In  this  chapter,  the  basic 
concepts  and  definitions  required  for  characterization  and  analysis  of  frequency  stability 
of  general  precision  oscillators  is  briefly  described. 


Clocks  are  the  instruments  we  use  to  measure  time  or,  to  be  more  precise, 
increments  of  time  elapsed  since  an  arbitrary  chosen  initial  reference.  Basically,  every 
clock  consists  of  two  different  parts:  an  oscillating  device  and  a  counter  that  accumulates 
the  number  of  oscillations  performed  since  it  is  initially  set  to  zero.  Thus,  the  initial 
reference  of  “time”  is  somewhat  arbitrary.  Clocks  have  evolved  from  astronomical  clocks 
based  on  planetary  movements  and  governed  by  the  laws  of  celestial  mechanics  to  clocks 
based  on  pendulums,  mechanical  springs,  or  the  piezoelectric  properties  of  quartz 
crystals.  Precision  atomic  clocks  used  today  are  based  on  the  quantum  transitions  of 
electrons  occurring  between  the  different  states  of  certain  elements  like  rubidium  or 
cesium  [1]. 


The  initial  definition  of  the  second  is  based  on  the  mean  solar  day,  which  has 
exactly  86400  seconds.  The  definition  of  the  mean  solar  day  is  based  on  astronomical 
computations.  As  new  precise  instruments  were  available  to  measure  increments  of  time, 
it  was  soon  noticed  that  there  were  irregularities  in  the  rotation  of  the  Earth;  in  particular, 
the  speed  of  rotation  was  observed  to  be  slowing  down  at  a  rate  in  the  order  of  one 
second  per  year.  It  was  not  until  1960  in  the  11th  General  Conference  on  Weights  and 
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Measures  (CGPM,  Conference  Generale  des  Poids  et  Mesures )  when  a  new  definition 
based  on  the  tropical  year  and  proposed  by  the  International  Astronomical  Union  was 
adopted.  Finally,  in  1967  in  the  13th  General  Conference  on  Weights  and  Measures  the 
definition  based  on  astronomical  theories  was  abandoned.  A  new  definition  based  on  the 
physical  properties  of  a  particular  isotope  of  cesium  was  adopted  and  confirmed  in  1997 
by  the  International  Committee  for  Weights  and  Measures  (CIPM,  Comite  International 
des  Poids  et  Mesures).  The  new  definition  was: 


The  second  is  the  duration  of  9  192  631  770  periods  of  the  radiation 
corresponding  to  the  transition  between  the  two  hyperfine  levels  of  the  ground 
state  of  the  cesium  133  atom. 


Even  atomic  clocks  based  on  cesium  present  differences  in  the  time  scale  that 
each  clock  generates.  Their  readings  are  different,  and  the  time  they  take  to  perfonn  a 
complete  oscillation  is  different  from  oscillation  to  oscillation  and  from  clock  to  clock. 
The  objectives  in  this  thesis  include  the  study  of  the  different  behavior  of  atomic  clocks 
in  order  to  assign  them  a  measure  of  performance  indicating  how  well  each  clock  is 
performing.  This  measure  of  performance  associated  with  each  clock  can  later  be  used  as 
a  weighting  factor  to  build  a  new  composite  time  scale  derived  from  the  data  supplied 
from  the  clocks  belonging  to  the  ensemble. 


These  weighting  factors  are  time  varying  numbers  that  contain  infonnation  of  the 
goodness  of  the  scale  generated  for  each  clock.  The  more  stable  the  frequency  generated 
by  a  clock,  the  more  weight  this  clock  will  have  in  the  composite  scale.  The  fact  that 
there  is  no  reference  to  compare  to  and  that  the  only  available  data  are  relative  measures 
among  real  clocks  of  similar  characteristics  make  this  problem  interesting  as  well  as 
challenging  [2]. 
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B. 


CLOCK  MODEL  AND  DEFINITIONS 


Every  atomic  clock  provides  as  outputs  different  electronic  signals,  some  of  them 
are  digital  signals  with  square  pulses,  while  some  are  continuous  like  the  sinusoidal 
signals  of  different  frequencies  (generally  1,  5  and  10  MHz)  that  can  be  represented 
according  to  the  following  model  [3,4]: 

F(O  =  F0(O-sin(®(O)  (1-1) 

Here  V(t)  represents  the  voltage  level  of  the  signal,  V0(t)  represents  the  magnitude  or 
amplitude  of  the  oscillations  and  <f>(/)  is  the  total  oscillator  phase,  starting  from  some 
arbitrary  origin  ®(t  =  0) .  We  will  assume  that  the  amplitude  V0(t)  for  atomic  clocks  can 
be  treated  as  a  constant  Vo. 


The  total  oscillator  phase  <f>(/)  can  be  decomposed  as  the  sum  of  a  purely 
deterministic  component  that  is  responsible  for  the  oscillation  at  a  fixed  nominal 
frequency  v0,  plus  a  stochastic  component  known  as  the  residual  phase  or  phase 
deviation  (ft): 


d>{t)  =  2-n  -v^-t  +  ^f)  (1.2) 

The  instantaneous  angular  frequency  co{t)  of  the  oscillations  is  defined  as  the 
time  derivative  of  the  total  oscillator  phase: 

®(0  =  n  ‘V0‘t  +  0(0)  =  2- K -H)+  fat)  (1.3) 

dt  dt 

Finally,  the  instantaneous  frequency  v(t)  is  defined  as: 


v(t)  = 


- — a(t)  =  v0  + 
2-n 


fat) 

2-n 


(1.4) 


This  quantity  can  be  considered  as  a  nominal  frequency  of  oscillation  v0  affected 
by  the  perturbation  term  <f>(t)/(2-n) .  This  is  the  quantity  we  are  interested  in,  since  it 
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defines  the  stability  of  a  clock.  The  instantaneous  frequency  of  the  oscillations  v(t)  may 
change  erratically  because  of  the  processes  involved  in  the  generation  of  the  signal  inside 
the  clock  and  also  because  of  external  factors  (environmental,  geographic,  etc.). 
Unfortunately,  since  infinite  bandwidth  measurements  are  not  possible,  an  exact  value  for 
the  instantaneous  frequency  is  not  measurable. 


The  statistical  quantity  that  is  used  in  stability  characterizations  of  oscillators  is  an 
adimensional  quantity  that  measures  the  relative  instantaneous  deviation  in  frequency 
from  its  nominal  or  predefined  value  v0 .  It  is  called  the  fractional  frequency  or  fractional 
frequency  fluctuation  and  is  defined  as 


v(t)  =  v^~v°  =  M 

V„  2 -7T-V0 


(1-5) 


For  the  same  reason  given  above  related  to  the  need  of  infinite  bandwidth,  it  is  not 
possible  to  measure  the  instantaneous  fractional  frequency,  but,  as  will  be  shown,  an 
average  value  can  be  directly  estimated  from  differences  in  time  between  clocks. 


The  time  deviation  x(t)  is  defined  as  the  integral  of  the  fractional  frequency  with 
respect  to  time.  It  has  dimensions  of  time  and  represents  the  difference  in  time  between 
two  clocks  or  the  time  difference  between  a  real  clock  and  an  ideal  one  oscillating  exactly 
with  a  frequency  of  v0  and  initially  synchronized  at  time  t  —  0 .  This  is  the  magnitude  that 
is  measured  in  time  laboratories  and  used  to  estimate  the  average  fractional  frequency. 


x{t)  =  \y(t)-dt 
0 


y(t )  = 

2-7t-v0  dt 


2  ■  n  ■  v0 


(1.6) 


Time  deviation  x(t)  and  phase  deviation  <f)(t)  are  linearly  related  and 
consequently  present  the  same  spectral  characteristics. 


4 


c. 


CLASSIFICATION  OF  NOISE  PROCESSES 


The  fractional  frequency  deviation  y(t) ,  the  time  deviation  x(t)  or  the  phase 
deviation  (j>(f)  are  all  affected  by  two  different  categories.  First,  there  is  a  systematic  or 
deterministic  component  that  is  generally  produced  by  two  situations  [5,6]: 

•  A  constant  frequency  of  oscillation  that  is  different  than  the  nominal 
value  vn . 

•  A  constant  frequency  drift  D,  when  the  frequency  of  a  clock  instead  of 
being  constant  changes  linearly  with  time. 


The  second  category  is  caused  by  random  processes  involved  inside  the  clocks. 
This  is  the  most  difficult  category  to  characterize  the  time  scale  generated  by  each  clock, 
since,  even  if  the  frequency  of  the  clock  is  changing,  as  long  as  it  changes  according  to  a 
detenninistic  or  predefined  pattern,  the  scale  can  always  be  compensated  by  algebraic 
manipulations.  On  the  other  hand,  the  random  effects  are  more  difficult  to  measure  and 
compensate.  Study  of  these  effects  constitutes  the  central  part  of  this  work. 


For  a  clock  presenting  an  initial  frequency  offset  and  a  constant  drift,  the  time 
deviation  can  be  expressed  as  [5]: 

x(t)  =  x0+y0-t  +  ^D-t2  +£(t)  (1.7) 

Here  x0  is  the  synchronization  error,  and  v0  is  the  fractional  frequency  offset  and  includes 
the  effects  produced  by  an  initial  frequency  of  oscillation  different  than  the  nominal  v0  at 

time  t  =  0 .  The  term  involving  the  drift  rate  D  includes  the  effects  of  a  linear  rate  of 
change  in  the  instantaneous  frequency  of  oscillation  v(t)  and  e{t)  represents  the  random 
effects.  In  Figure  1,  the  effect  of  frequency  offset  and  frequency  drift  in  the  quantities 
x(t)  and  y(t)  is  shown.  In  the  figure,  the  two  systematic  error  cases  are  presented:  a 


5 


clock  with  a  constant  frequency  but  different  than  the  nominal  v0  and  another  clock  with 

a  negative  frequency  drift.  The  effects  produced  by  environmental  fluctuations,  humidity, 
temperature,  shock,  vibrations,  radiation  etc.,  can  also  induce  other  types  of  systematic 
deviations  (i.e.  modulation  sidebands)  that  will  not  be  treated  here  [5,6]. 


Fractional  Frequency  Deviation 


Time  Deviation 


Figure  1.  Deterministic  frequency  and  time  deviations 


The  types  of  noise  present  in  atomic  clocks  are  classified  according  to  their 
spectral  characteristics.  In  the  frequency  domain,  power-law  spectral  densities  serve  as  a 
reasonable  and  accurate  model  of  the  random  fluctuations  in  precision  oscillators.  In 
practice,  these  random  fluctuations  are  often  represented  by  the  sum  of  five  such  noise 
processes  assumed  to  be  independent,  as  [3,7,8]: 


Syif) 


'  a=2 

0  </</* 

a=- 2 

0 

f>fh 

(1.8) 
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Where  S  (f)  is  the  power  spectral  density  (PSD)  of  the  fractional  frequency, 
ha  is  a  constant  associated  with  each  noise  type  and  fh  is  the  high  frequency  cut  off 
frequency.  Sy  (/)  is  proportional  to  the  squared  magnitude  of  the  Fourier  Transform  of 
the  fractional  frequency. 


The  duality  time-frequency  and  the  relationships  in  the  time  domain  among  the 
fractional  frequency  y(t) ,  time  deviation  x(t)  or  phase  deviation  (f>{t)  have  a  counterpart 
in  the  frequency  domain,  so  it  is  frequent  to  find  in  the  literature  equivalent  expressions 
for  the  noise  models  as  a  function  of  the  power  spectral  densities  of  the  time  deviation 
Sx  (/)  or  the  phase  deviation  S,  (/) .  All  of  them  are  equivalent.  The  dual  relationships 

of  (1.6)  in  the  frequency  domain  are  [9]: 

w)=7 — s,(/)=4w)  <L9> 

(2-n-v„Y  v0 

The  noise  processes  present  are  usually  identified  by  plotting  in  logarithmic  axes 
the  PSD  and  then  fitting  a  first  degree  polynomial.  The  corresponding  slope  of  the 
straight  line  relates  it  to  the  noise  type. 


Table  1  summarizes  five  commonly  used  models  for  noise  and  their 
characterizations  in  the  frequency  domain  [3,6,10]: 
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Frequency 

Noise  Process 

Domain 

sy(f) 

Sx(f)  or  S,(f) 

Slope:  log(PS'D)  vs.  log (/)  a,  /? 

a 

P 

Random  Walk  Frequency  Modulation 

-2 

-4 

Flicker  Frequency  Modulation 

-1 

-3 

White  Frequency  Modulation 

0 

-2 

Flicker  Phase  Modulation 

1 

-1 

White  Phase  Modulation 

2 

0 

Table  1 .  Noise  processes  used  in  modeling  frequency  instability  (frequency) 

D.  TIME  DOMAIN  CHARACTERIZATION,  THE  ALLAN  VARIANCE 

1.  Definition  and  Estimate. 


The  quantity  of  interest  is  the  fractional  frequency,  or  to  be  more  precise,  the 
residual  fractional  frequency  which  corresponds  to  the  same  concept  once  the 
deterministic  components  or  systematic  effects  have  been  removed  [10,1 1]. 


The  fractional  frequency  data  provides  information  about  the  instantaneous 
frequency  of  the  oscillator.  If  it  has  a  non  zero  but  nearly  constant  value,  it  means  that  the 
oscillator  is  generating  a  nearly  constant  (stable)  signal  but  at  a  frequency  different  than 
the  nominal  v0  at  which  it  is  supposed  to  oscillate.  This  is  not  a  bad  sign  as  long  as  it  has 

little  variation.  The  parameter  that  characterizes  the  stability  is  consequently  a  function  of 
the  variability  of  the  fractional  frequency  y(t)  . 
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The  instantaneous  value  of  the  fractional  frequency  of  one  clock  compared  to 
another  cannot  be  computed,  but  an  average  value  during  a  time  interval  can  be  obtained 
from  the  time  differences  between  these  clocks.  Let  us  suppose  that  x( t )  represents  the 
time  difference  between  two  clocks  at  time  t.  We  usually  compare  the  clocks  at  equally 
spaced  instants  of  time,  that  is,  the  values  of  x{t)  known  to  us  are  xm  =  x(m  ■  r0)  where  m 

takes  integer  values  and  r0  represents  the  time  spacing  or  sample  period  between 
measurements.  The  time  series  containing  the  measurement  of  time  differences  is 
represented  by  |x0,  jq,  x2,  x3, ...,  jc^}  ,  where  N  is  the  number  of  time  difference  data 
points  available. 


The  relationship  between  the  fractional  frequency  and  the  time  difference  is 
introduced  in  (1.6);  average  values  of  fractional  frequency  y(t)  can  be  estimated  by 
approximating  the  derivative  with  the  ratio  of  increments: 


y(t) 

m 


dx{t ) 
dt 

A x{t) 
At 


x{t  +  t)  -  x{t) 

T 


(1.10) 


During  periods  of  time  of  r0  seconds,  average  values  for  y(t)  can  be  computed: 


v(0-  r0  <  t  <  1  •  r0)  =  y0  = 


y(l-r0  <t<2-r0)  =  yl  = 


x(\-t0)-x(Q-t0) 

L) 

x(2-t0)-x(1  -r0) 


y((N  -  2)  •  t0  <  t  <  (N  - 1)  •  r0)  =  yM_\ 


x((N  - 1)  •  r0)  -  x((N  -  2)  •  r0 )  _  xN_x  -  xN_2 
T>  To 

(1.11) 


Here  M  =  N  - 1  is  the  number  of  fractional  frequency  data  points.  The  time  series 
obtained  with  the  averaged  fractional  frequencies  is  represented  by 
y{t)  =  {y0,yl,y2,yi,y4,...yM_l),  and  will  provide  the  required  data  to  characterize  the 

stability  associated  with  a  particular  pair  of  clocks.  This  is  a  very  straightforward 
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procedure  since  the  time  deviations  x(m-  r0)  are  just  the  differences  in  time  between 
clocks,  which  can  be  measured  with  a  counter  and  with  negligible  measurement  noise  (on 
the  order  of  picoseconds). 


If  M  values  of  averaged  fractional  frequency  y(  are  computed  using  (1.11),  there 
are  many  ways  to  analyze  these  data.  Historically,  people  have  used  the  sample  standard 
deviation  as  a  measure  of  variability  for  the  statistic  yi  and  estimated  by: 

I  i 

°std.  dev  =  y  y  •  £  ( yt  -  y  f  ( 1  • 1 2) 

The  usual  convention  of  using  the  same  symbol  for  the  definition  and  for  the 
estimator  of  quantities  is  followed  in  this  thesis,  except  that  estimators,  as  in  expression 
(1.12),  will  be  represented  with  a  “hat”  to  avoid  confusion.  In  the  last  expression,  y 
represents  the  average  fractional  frequency  over  all  the  data  set. 


This  approach  presents  several  difficulties;  in  particular,  depending  on  the 
spectral  components  present  in  the  noise,  the  classical  or  sample  standard  deviation  may 
not  converge  to  a  fixed  value  [11],  (it  only  converges  for  White  Noise  Frequency 
Modulation  WFM).  In  general,  it  is  a  function  of  the  number  of  time  data  points  N,  which 
is  not  acceptable.  For  example,  for  two  clocks  that  have  slightly  different  stable 
frequencies  of  oscillation,  the  standard  deviation  of  the  fractional  frequency  fluctuation  is 
a  monotonically  increasing  function  of  the  sample  size  and  increases  without  limit  (see 
Figure  2).  It  is  clear  that  some  other  measure  of  variability  of  the  fractional  frequency 
must  be  chosen.  The  IEEE  has  adopted  a  standard  measure  known  as  the  “Allan 
variance”[3]  also  known  as  two-sample  variance. 


Figure  2,  obtained  from  [10],  shows  the  ratio  of  the  usual  sample  variance  to  the 
Allan  variance  as  a  function  of  the  number  of  samples  N  for  the  five  classes  of  noise 
present  in  clocks.  It  shows  why  the  sample  variance  is  not  a  convenient  measure  of 
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frequency  stability;  i.e.  the  number  of  samples  in  the  computation  has  a  very  significant 
effect  in  the  final  value  of  the  estimate,  the  lack  of  convergence  is  also  evident  except  for 
the  White  Noise  Frequency  Modulation. 


Figure  2.  Ratio  of  sample  /Allan  variance  for  different  noise  processes  (After  ref.  [10]) 

The  square  root  of  the  Allan  variance  is  called  Allan  deviation  and  is  defined  by: 

/i  a1/2 

<Ty(T)  =  ( -\y(t  +  T)-y(i)f\  ,T>  0  (1.13) 

Where  the  symbols  “ /  )  ”  denote  infinite  time  average  and  the  function  y(t)  is 
represented  by  the  piecewise  linear  function  obtained  from  the  series 
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{ v, ,  y2,  y3,  y4,  ...yM}  computed  at  intervals  of  z  time  units  as  in  (1.1 1).  The  factor  of  Vi 

appearing  in  the  definition  (1.13)  is  introduced  only  for  normalization  purposes,  when 
there  are  only  two  data  samples,  the  classical  standard  deviation  and  the  Allan  deviation 
have  exactly  the  same  value  as  shown  in  Figure  2,  no  matter  what  kind  of  noise  process  is 
present.  In  practice,  crv(r0)  can  be  estimated  from  a  finite  set  of  M  fractional  frequency 

data  points  as  follows: 


<7y(h) 


1 

P-(M-I) 


M- 1 

Z(k,+i-y,)2 


i= 1 


(1.14) 


Again,  the  estimator  has  been  represented  with  the  “hat”  notation.  The  relationship 
between  the  Allan  deviation  and  the  classical  Standard  deviation  as  a  function  of  the 
sample  size  N  for  the  common  types  of  noise  is  shown  in  Table  2: 


Noise  Type 


Standard  deviation  cr 


std.  dev 


Random  Walk  Frequency  Modulation 


cr 


std.  dev  y 


(  r  iN 


Flicker  Frequency  Modulation 


® std .  dev  ®y  (7)  ' 


N  ■  log(A) 

2  •  (A  - 1)  •  log(2) 


White  Frequency  Modulation 


&  std.  dev  (7) 


Flicker  Phase  Modulation 


° std.  dev  *  ®y  (7)  ' 


2-(A-l) 
1  3  •  N 


White  Phase  Modulation 


° std.  dev  17 y  (7)  ' 


2-(N-l) 
3  •  N 


Table  2.  Relationship  between  classical  standard  deviation  and  Allan  deviation  (After  ref. 

[18]) 


The  same  concept  of  Allan  variance  can  be  generalized  to  include  Allan 
covariance  between  two  time  series.  For  two  given  time  series  of  fractional  frequencies 


12 


y\t)  and  y2(t )  defined  according  to  (1.11),  the  Allan  covariance  of  the  two  series  is 
defined  by: 

(r)  =  ■  (y  {t  +  T)-yx{t))\y2{t  +  z)-y2{t))^,T>0  (1.1 5) 

which  can  also  be  estimated  by: 

i  M- 1 

*yfo>)  =  2-(M  -  (L16) 


There  is  an  alternative  definition  of  the  Allan  deviation  and  Allan  covariance  as  a 
function  of  the  time  deviation  series,  x( I ) ,  instead  of  the  fractional  frequency.  For  the 
Allan  deviation  it  can  be  derived  from  (1.13)  and  (1.10): 


1/2 


o-y(t)  =  < 


,2-r 


:{-x(t  +  2-T)  +  2-x(t  +  r)-x(t))j  ,  r>0  (1.17) 


The  corresponding  estimator  is: 


N-2 


°y(? o)  =  . 


'2-(A-2)-r( 


Z(-x;+2+2  -xM-Xi) 


(1.18) 


0  1=1 


Expressions  (1.14)  or  (1.18)  estimate  the  Allan  variance  for  an  integrating  time 
equal  to  the  sampling  period  r0 ,  but  values  of  the  Allan  variance  for  integrating  times 

multiples  of  r0  are  also  needed  because,  as  will  be  shown  later,  they  allow  us  to  estimate 
the  type  of  noise  process  present  in  the  data.  The  process  for  computing  values  of 
<j2v(t  =  m-T0)  at  multiples  m  of  the  fundamental  sampling  period  r0  is  straightforward. 

Starting  from  the  initial  time  series  containing  the  fractional  frequency  values 
{?,,  y2,y3,  y4,...yM}  ,  to  compute  the  Allan  variance  or  deviation  for  an  integrating  time 

of  t  =  m  ■  r0  it  is  only  required  to  obtain  a  new  time  series  derived  from  the  original 
where  each  tenn  is  an  average  of  a  subset  of  m  terms  of  the  original  data  set  and  is 
computed  as  [1 1]: 


» 


y m-i  y m-i+l  T m-i+2  *"  T m-i+m— 1 


m 


(1.19) 
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Once  the  new  time  series  {y0 yt y2 y3 y4 }  has  been  obtained,  the 
expression  given  in  (1.14)  can  be  used  to  find  &v(m  ■  t0)  ,  in  this  case,  M'  represents  the 
number  of  tenns  in  the  new  series.  This  process  can  be  done  for  any  multiple  of  r0 ,  but 
in  practice,  this  computation  is  usually  done  for  multiples  that  are  powers  of  2,  i.e. 
m  =  1,  2,  4,  8, ... . 


Because  of  the  averaging  process  involved  in  the  computation  of  the  new  time 
series  when  m  ^  1 ,  as  the  number  m  increases,  there  are  less  data  points  to  compute 
av(ni  ■  r0)  and  consequently  the  confidence  interval  for  cr  (in  ■  r0)  also  increases.  To 

minimize  this  negative  effect,  modern  algorithms  make  a  better  use  of  the  data  available 
obtaining  smaller  confidence  intervals.  The  standard  algorithm  used  to  compute  the 
statistic  is  the  overlapped  algorithm  and  uses  the  time  deviation  series  or  raw  data 
obtained  from  the  clock  comparison:  jc.,  Vz  =  {0,1,2,..., N-l)  ,  the  overlapped  estimate  is 
given  by  [3]: 

Ji  N—2-m 

?  7T7  9 - N - 2~2  X  Z  (  Xi+2-m  +  2  •  Xi+m  -  x,  )2  (  1 .20) 

2  ■  [N  -  2  ■  m)  ■  m  -r0 

The  fact  that  the  Allan  variance  or  Allan  deviation  does  converge  for  the  five 
types  of  noise  present  in  precision  oscillators  makes  it  suitable  to  measure  the  variability 
of  the  fractional  frequency.  However,  it  is  not  the  only  statistic  that  shows  this  property, 
other  variances  widely  used  with  similar  properties  include  The  Modified  Allan  variance 
[12,13],  The  Total  variance[l  1],  The  Time  variance[14],  The  Hadamard  variance  (or 
three-sample  variance)[15]  and  others  [16,17]. 


The  Allan  variance  also  presents  some  inconveniences;  the  most  important  is  that 
by  using  the  Allan  variance  is  not  possible  to  differentiate  White  Phase  Modulation 
(WPM)  noise  from  Flicker  Phase  Modulation  (FPM)  noise.  When  those  types  of  noise  are 
present  and  a  more  accurate  characterization  is  needed,  then  other  statistics  that  are 
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capable  of  discriminating  these  two  processes  can  be  used,  like  the  modified  Allan 
variance  [12]. 

2.  A  Simple  Example  to  Compute  Allan  Variance. 


To  illustrate  the  computation  of  the  Allan  variance  or  Allan  deviation,  a  simple 
example  is  provided  below.  Let  us  suppose  that  the  time  differences  between  two  clocks 
have  been  measured  every  second  and  nine  measurements  have  been  made.  The  real  time 
difference  between  them  x(t)  is  a  continuous  function  represented  in  Figure  3,  but  only 
nine  data  samples  are  known. 


Time  difference  between  two  clocks 


Figure  3.  Sample  data  to  compute  Allan  variance 
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The  interval  between  samples  is  r0  =  1  second,  so  the  Allan  variance  can  be 
computed  at  multiples  of  this  r0 :  r  =  n  •  r0,  n  =  {1,2,3,...} .  The  computation  for  r  =  Is 
is  given  in  Table  3. 


K 

xk 

(xl(T6) 

W  =(At+ !-**)/* 
(xl(T6) 

yk+x-yk 

(xio^6) 

(yk+ 1  -ykf 

(xlO12) 

°y(?  =  1) 

(xlO6) 

0 

x0  =0 

To  =43. 6 

2.5 

6.25 

1 

Xj  =  4  3 . 6 

Ti  =46.1 

-14.2 

201 . 64 

2 

x2  =  8  9 . 7 

v2  =  3 1 . 9 

10.2 

104 . 04 

3 

x3  =  1 2 1 . 6 

T3=42.1 

2 . 6 

6.76 

4 

x4  =  1 6  3 . 7 

r- 

II 

l?s 

-5.1 

26.01 

5.67 

5 

x5  =  2  0  8 . 4 

T5  =  3  9  •  6 

1.4 

1 . 96 

6 

x6  =  2  4  8 . 0 

T6=41.0 

-10.2 

104 . 04 

7 

x7  =  2  8  9 . 0 

y7  =30.8 

— 

— 

8 

x8  =  3 1 9 . 8 

— 

— 

— 

Table  3.  Steps  for  computing  Allan  deviation  for  v  =  Is 


Where  <tv(t  =  1)  is  computed  as: 

<r,(r  =  l  )  =  ]^i(y,.,-y,)2  =5.67-10  -6s  (1.21) 

V  2  '  /  i  =0 

For  values  of  r  different  than  r0 ,  the  computation  follows  a  similar  process. 
However,  the  classical  estimate  does  not  use  all  of  the  samples  of  time  differences 
available.  The  steps  to  compute  av(r  =  2)  is  given  in  Table  4: 
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K 

(xl(T6) 

yk  =i.xk+2-xk)lr 
(xl(T6) 

yk+2-yk 

(xio-6) 

(yk+2-yk)2 

(xlCT12) 

ay  (T  =  2) 

(xlO6) 

0 

x0  =  0 

ll 

4^ 

4^ 

00 

on 

-7 . 85 

61 . 62 

1 

Xj  =  4  3 . 6 

— 

— 

— 

2 

x0  =  8  9 . 7 

y2  =37.00 

5.15 

26.52 

3 

*3=121.6 

— 

— 

— 

4 

x4  =  1 6  3 . 7 

V4=  42.15 

-6.25 

39.06 

4.60 

5 

x5  =  2  0  8 . 4 

— 

— 

— 

6 

x6  =  2  4  8 . 0 

v6=35.90 

— 

— 

7 

x7  =  2  8  9 . 0 

— 

— 

— 

8 

xg  =  3 1 9 . 8 

— 

— 

— 

Table  4.  Steps  for  computing  Allan  deviation  for  z  =  2s 


Where  a  (z  =  2)  is  computed  as: 

^  =  2)  =  jA£(y!M_y!t);  =4.60- 10^  (1.22) 

A  better  algorithm  to  estimate  the  Allan  variance  for  values  of  z  different  than  z(] 
should  use  all  available  fractional  frequency  data.  This  is  what  the  overlapped  estimate 
does.  It  is  the  algorithm  used  in  this  thesis.  Both  estimators  are  unbiased  for  the  Allan 
variance,  but  in  the  case  of  the  overlapped  version,  the  confidence  interval  is  smaller. 
Confidence  intervals  for  the  Allan  variance  will  be  treated  in  the  next  section.  The 
overlapped  algorithm  to  compute  the  Allan  variance  is  shown  in  Table  5: 
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Table  5.  Steps  for  computing  (overlapped)  Allan  deviation  for  r  =  2s 


Where  cry(z  =  2)  is  computed  as: 

°-,(r  =  2  )  =  ,fl7E(yM-yi)2  =3.95-10“*  (1.23) 

\  k=  0 

3.  Confidence  Intervals  for  the  Allan  Variance. 

The  simplest  method  to  compute  the  confidence  intervals  for  cr  (r)  [3,6,1 1], 

which  assumes  a  symmetric  Gaussian  distribution,  uses  the  following  expression  to 
compute  95%  confidence  intervals  for  the  Allan  variance  [19]: 

I  a  =  1-96  •  cr„(r)  •  Ka  /  y[N  (1.24) 
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where  Ia  is  the  half  amplitude  of  the  interval,  Ka  is  a  constant  related  to  the  type  of  noise 

present  and  N  is  the  number  of  data  points  used  in  the  estimate  and  a  corresponds  to  the 
slope  in  logarithm  axes  of  log  (PSD)  vs.  log(/)  . 


The  values  for  Ka  for  the  different  noise  types  are  given  in  Table  6. 


Noise  Type 

RWFN 

FFM 

WFM 

FFP 

WPM 

a 

-2 

-1 

0 

+  1 

+2 

Ka 

0.75 

0.75 

0.87 

0.99 

0.99 

Table  6.  Values  of  Ka  to  compute  confidence  intervals  for  <j2 


Another  approach  to  computing  confidence  intervals  is  based  on  the  distribution 
of  the  estimate  for  the  Allan  variance.  The  Allan  variance  estimator  is  distributed 
according  to  a  Chi-Squared  distributional  1]: 


2 


d.f.xs2 


(1.25) 


Where  s2  is  the  sample  Allan  variance,  d.f.  corresponds  to  the  degrees  of  freedom 
and  cr2  is  the  true  value  of  the  Allan  variance.  The  unknown  parameter  in  the 
computation  is  the  degrees  of  freedom  of  the  distribution.  Much  research  has  been  done 
in  the  past  regarding  the  number  of  degrees  of  freedom.  Both  analytical  and  empirical 
expressions  have  been  proposed[6].  The  number  of  degrees  of  freedom  depends  on  the 
algorithm  used  to  compute  the  Allan  variance  and  more  importantly  on  the  type  of  noise 
process  present  in  the  clocks.  In  this  thesis,  the  algorithm  used  to  compute  the  Allan 
variance  is  the  overlapped  version  since  it  provides  smaller  confidence  intervals  for  the 
same  number  of  data  points. 


The  empirical  expressions  for  the  number  of  degrees  of  freedom  that  are 
recommended  in  characterization  of  precision  oscillators  are  given  in  Table  7  for  the 
overlapped  Allan  variance. 
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Noise  type 

Degrees  of  Freedom 

Random  Walk  Frequency  Modulation 

(RWFM) 

N —  2  (iV  — 1)~  —3  •  m  •  ( A  —  l)  +  4  •  m2 
m  ''  (A-3)2 

Flicker  Frequency  Modulation 

(FFM) 

2-(A-2)2 

V  ;  (m  =  1) 

2.3 -TV -4.9 

5-N2 

,  ,  (m>  2) 

4-m-(N  +  3  ■  w) 

White  Frequency  Modulation 

(WFM) 

4 -m2  f  3-(A-l)  2-(A-2)^| 

4-m2  +5  ^  2-m  N  J 

Flicker  Phase  Modulation 

(FPM) 

1  f,  fN-l']  ,  ((2-m  +  l)x(N-l)\\ 

exp  log  xlog  „ 

V  l  V2‘mJ  V  4  )) 

White  Phase  Modulation 

(WPM) 

{N +  \)x{N -2-m) 

2  -{N  —  m) 

Table  7.  Empirical  expressions  for  the  degrees  of  freedom  of  overlapped  Allan  variance 

(After  ref.  [3]) 


4.  Time  -  Frequency  Duality. 


The  characterization  of  precision  clocks  or  oscillators  can  be  done  in  both 
domains:  time  or  frequency.  The  noise  processes  that  may  be  present  are  defined  in  the 
frequency  domain  because  this  domain  provides  an  adequate  and  straightforward 
definition  by  using  the  Power  Spectral  Density  function.  However,  the  duality  between 
both  domains  permits  us  to  identify  the  type  of  noise  present  while  working  exclusively 
in  the  time  domain.  That  is,  the  plot  of  the  Allan  variance  or  Allan  deviation  for  different 
values  of  r  almost  completely  characterizes  the  process  [3]  (for  a  complete 
characterization  in  the  time  domain  the  modified  Allan  variance  could  be  used  instead). 
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When  the  Allan  variance  is  plotted  on  logarithm  axes,  the  slope  of  the  curve  //  is 


related  to  the  noise  type  according  to  Table  8  [3,6,10]: 


Noise  Process 

Time 

Domain 

S' 

Mod  cr 

y 

Slope:  log(crv(r)) vs. log(r)  ju,ju' 

A 

fill 

A’ 

Random  Walk  Frequency  Modulation 

1 

1/2 

1/2 

Flicker  Frequency  Modulation 

0 

0 

0 

White  Frequency  Modulation 

-1 

-1/2 

-1/2 

Flicker  Phase  Modulation 

-2 

-1 

-1 

White  Phase  Modulation 

-2 

-1 

-1/2 

Table  8.  Noise  processes  used  in  modeling  frequency  instability  (Time) 


Complete  characterization  in  the  time  domain  is  not  possible  using  only  the  Allan 
variance.  In  Table  8,  White  Phase  Modulation  and  Flicker  Phase  Modulation  show  the 
same  slope.  However,  the  modified  Allan  variance  provides  a  way  to  separate  them,  as 
shown  in  the  fourth  column. 


21 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


22 


II.  CHARACTERIZATION  OF  THE  CLOCK  ENSEMBLE 


A.  STATEMENT  OF  THE  PROBLEM 


Each  clock  has  its  own  intrinsic  characteristics  that  have  an  immediate  effect  on 
the  stability  of  the  time  scale  that  it  generates.  Unfortunately,  it  is  not  possible  to  perform 
absolute  measurements  on  a  particular  clock  because  there  is  not  a  perfect  clock  to 
compare  to.  We  can  only  measure  the  time  phase  difference  between  pairs  of  real  clocks. 
In  other  words,  we  only  have  access  to  relative  measurements  that  involve  the 
instabilities  generated  inside  both  clocks. 


The  main  effort  in  this  chapter  is  to  find  an  absolute  measure  of  performance 
associated  with  each  clock  in  order  to  assign  a  weight  that  will  later  be  used  to  build  an 
ensemble  scale.  The  number  of  clocks  belonging  to  the  ensemble  is  represented  by  K. 
The  available  data  to  conduct  this  study  is  data  obtained  in  laboratory  consisting  of  K- 1 
time  differences  measured  between  all  K  clocks  in  the  ensemble  with  respect  to  one 
arbitrary  clock  that  is  chosen  as  a  reference.  These  time  differences  are  measured  at 
equally  spaced  instants  of  time  and  therefore  each  pair  of  clocks  generates  a  time  series. 


Here  is  an  important  observation  related  to  the  concepts  of  stability  and  its  change 
with  time.  In  the  rest  of  this  thesis,  the  quantity  “time”  is  simultaneously  playing  the  role 
of  both  dependent  and  independent  variable.  The  main  objective  of  this  work  is  devoted 
to  characterization  of  the  goodness  of  the  time  scale  generated  by  each  source  of  time  or 
clock.  The  assumption  that  there  is  not  a  perfect  time  scale  generator  seems  to  contradict 
the  idea  that  comparisons  between  real  clocks  are  taken  at  “equally-spaced”  instants  of 
time.  Obviously,  there  are  not  such  “equally-spaced”  instants  of  time  since  there  is 
instability  always  present,  but  it  is  important  to  understand  the  different  orders  of 
magnitude  that  are  related  to  “time”  when  it  is  playing  the  role  of  dependent  or 
independent  variable.  Typical  instabilities  of  atomic  clocks  are  in  the  order  of  a  few 
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nanoseconds  per  day,  the  nanoseconds  measure  the  instability  (dependent  variable)  while 
the  day  is  the  order  of  the  independent  variable  and  it  makes  no  significant  difference  if 
the  length  of  the  day  is  24  hours  or  24  hours  ±  a  few  nanoseconds. 


B.  CONSTRUCTING  A  COMPOSITE  CLOCK 


One  of  the  final  goals  in  the  characterization  of  the  ensemble  of  clocks  is  the 
construction  of  an  optimum  time  scale  with  respect  to  the  stability  in  the  frequency  that  it 
generates.  In  order  to  do  that,  the  time  varying  behavior  of  each  clock  must  be  analyzed 
and  an  ensemble  weight  for  each  clock  must  be  computed. 

Let  y  =(vvy2,...,yK)T  be  a  vector  of  averaged  fractional  frequencies 

corresponding  to  an  averaging  time  t  for  the  K  clocks  belonging  to  the  ensemble.  Let 
z  =  wy  be  a  scalar  composite  fractional  frequency  where  w  is  the  vector  of  weights 

given  by  w  =  ( w, ,  w2, ...,  wK )  .If  y  is  multivariate  normal  with  mean  0  and  covariance 
matrix  if,  what  should  the  weights  be  in  order  to  minimize  the  variance  of  z? 


The  traditional  way  to  assign  a  weight  to  each  clock  is  based  on  the  variance 
tenns  of  if,  that  is,  the  elements  in  the  main  diagonal.  Since  if  represents  the  covariance 
matrix  of  the  fractional  frequency  deviations,  the  bigger  the  variance  term  corresponding 
to  a  particular  clock,  the  more  variability  in  the  generated  frequency  and  consequently, 
the  less  weight  should  be  assigned  to  that  clock. 


The  mathematical  function  used  to  assign  weights  generates  them  proportionally 
to  the  inverse  of  the  corresponding  variance  term  subjected  to  the  condition  that  all  of  the 
weights  must  add  up  to  1.  For  example,  for  an  ensemble  of  three  clocks  with  covariance 
tenns  given  by  rn ,  r22  and  r33 ,  the  corresponding  weights  are  computed  by: 
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(2.1) 


l/rn+l/r22+l/r33 

1/^22 

W9  = - — - 

l/rn  +  l/r22  +  l/r33 

!/  r33 

w3  = - — - 

l/rn +1//-22 +1/7-33 


In  general,  for  an  ensemble  of  K  clocks,  the  weight  of  clock  i  can  be  obtained  by: 


j= 1 


Vie  {1,2,...,AT} 


(2.2) 


This  classical  approach  is  derived  from  the  Three  Cornered  Hat  method,  in  which  the 
clocks  are  assumed  to  be  uncorrelated  and  R  is  therefore  diagonal. 


When  correlations  between  clocks  occur,  R  is  no  longer  diagonal  and  a  more 
accurate  technique  can  be  used  instead.  If  the  weights  of  the  clocks  are  grouped  in  the 
vector  w ,  the  weights  can  be  obtained  in  general  by  minimizing  the  variance  of  z 
subjected  to  the  constraint  that  the  weights  must  add  up  to  1  for  the  solution  to  be 
unbiased: 


Min  Var(z)  =  w-R  ■  wT 

Wj 

s.t.  Bwt=  1 


(2.3) 


Here  the  vector  B  is  a  vector  of  K  ones:  B  =  (1, 1, 1, ...,  1) . 


Since  R  is  symmetric  and  non-singular,  it  can  be  shown  that  there  is  only  one 
stationary  point  for  the  function  w  ■  R  ■  wT  subject  to  B  ■  \vT  =  1  that  is  given  by: 

w*  =R  l  BT  (B  R1  B7)1  (2.4) 
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This  last  equation  can  be  used  whenever  R  is  no  longer  diagonal.  When  R  is 
diagonal,  both  equations  (2.2)  and  (2.4)  provide  the  same  results.  Since  the  correlation 
between  clocks  is  small  in  practice,  there  is  little  difference  between  the  two  approaches. 


In  any  case,  it  is  clear  that  knowledge  of  the  matrix  R  is  essential  to  the  problem 
of  constructing  a  composite  time  scale.  We  now  turn  to  the  problem  of  estimating  R. 


C.  MATHEMATICAL  FORMULATION 


Assume  for  the  moment  that  there  exists  an  ideal  clock  that  provides  a  stable  time 
scale  and  let  x'  be  a  stochastic  process  associated  with  the  7th  clock  of  the  ensemble.  Let 
now  x'm  be  a  realization  of  this  process  at  the  instance  of  time  tm .  The  interpretation  of  the 

stochastic  processes  x  can  be  thought  as  the  relative  time  deviation  between  the  7th  clock 
and  the  ideal  clock. 


Let  now  the  vector  x‘  =  \x[ :  k  =  1,2,...,  a|  be  the  sequence  of  N consecutive  and 

equally  spaced  realizations  of  the  process  x  .  This  vector  contains  the  relative  time 
deviation  (differences)  between  clock  i  and  an  ideal  clock  measured  at  N  different  equally 
spaced  instants  of  time. 


X  ^  Xj  ,  X2  ,  Xj  ,  .  .  .  ,  Xjy  ^ 


(2.5) 


The  expected  value  of  the  process  x'  can  be  estimated  as  a  function  of  its 
elements  as: 


x'  = 


—  —  (x[  +  x  2  +....  +  x'N )  —  —  'V' 
N{  1  2  Nj  Ntt 


(2.6) 


and  arranged  into  a  constant  vector  of  N  elements  as: 


x1  =  [x',xl,  xl, ...,  x‘ ) 


(2.7) 
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For  purposes  of  compactness,  the  expected  values  associated  with  all  of  the  K 
processes  for  the  clocks  in  the  ensemble,  can  also  be  arranged  in  a  NxK  matrix: 

X=(x\x2,  x3, ...,**)  (2.8) 

Likewise,  the  N  samples  or  time  differences  between  the  K  clocks  and  the  ideal 
one  can  also  be  arranged  into  a  similar  NxK  matrix: 

*  =  (*',  x2,  x\...,xK)  (2.9) 

If  the  K  processes  of  time  deviations  are  arranged  in  the  vector  u 
u  =  (x1,  x2,  x3, ...,  xk )  ,  then  the  covariance  matrix  associated  with  the  processes  can  be 
defined  as: 

R  =  e[u-ut}-E{u}-E[ut)  (2.10) 

where  the  operator  E  {  }  stands  for  mathematical  expectation.  Now,  the  covariance 

matrix  of  the  processes  associated  with  the  available  samples  of  time  deviations  for  the  K 
clocks  in  the  ensemble  can  be  estimated  by  [20,21,23]: 

•(*-*)  <21» 

The  KxK  matrix  R  is  referred  to  as  the  absolute  covariance  matrix  since  it  is  the 
covariance  matrix  involving  the  processes  relating  the  time  differences  between  real 
clocks  in  the  ensemble  with  respect  to  fictitious  ideal  one.  A  generic  element  of  the 
matrix  R  can  be  estimated  by: 

rtj  =^I-|(V  -  x1  f  -(xJ-  xJ )  (2. 12) 

When  i  =  j ,  r..  represents  the  estimate  of  the  variance  of  the  time  deviations  of 

clock  i  and  it  is  a  nonnegative  quantity;  when  /  ^  j ,  represents  the  estimate  of 

covariance  between  clocks  i  and  j  and  can  be  positive  or  negative.  The  accuracy  in  the 
estimation  of  R  improves  with  the  number  of  samples  available. 
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The  final  goal  here  is  to  estimate  the  absolute  covariance  matrix  R  in  order  to 
characterize  the  ensemble  of  clocks  and  to  build  a  composite  time  scale  using  all  of  the 
clocks  belonging  to  the  ensemble.  Unfortunately,  there  is  no  such  ideal  clock  that  has 
been  assumed  in  the  definition  of  the  above  quantities,  consequently  the  elements  r..  of 

matrix  R  can  not  be  directly  computed  or  estimated  because  it  is  not  possible  to  know  the 
time  differences  x'k  between  the  j  clock  and  the  ideal  one  at  a  particular  time  k. 


Let  us  suppose  that  the  ensemble  of  time  standards  consisting  of  K  clocks  of 
similar  characteristics  are  arranged  in  such  a  way  that  it  is  possible  to  measure 
periodically  the  time  differences  of  readings  between  each  clock  belonging  to  the 
ensemble  and  one  arbitrary  clock,  also  included  in  the  ensemble,  that  has  been  chosen  as 
a  reference.  Without  loss  of  generality,  let  us  assume  that  the  Xth  is  chosen  as  the 
reference.  At  each  particular  instant  of  time,  it  is  possible  to  measure  all  pairwise 
differences  between  the  K- 1  first  clocks  and  the  reference.  In  a  strict  sense,  it  is  not 
possible  to  perform  simultaneous  measurements  of  all  different  pairs  of  clocks.  However, 
since  the  time  required  to  measure  all  pair-wise  comparisons  among  all  clocks  is  in  the 
order  of  a  few  seconds  and  the  time  between  successive  measurements  is  in  the  order  of 
days,  the  effect  is  always  ignored  (this  can  produce  negative  effects  when  not  properly 
taken  into  account,  as  will  be  explained  below). 


It  is  then  possible  to  express  the  measured  data  as  a  function  of  the  differences  of 
time  between  each  clock  and  the  ideal  one,  that  is,  let  yf  be  the  time  difference  of 
readings  between  clocks  i  and  j  at  time  k,  then: 

yf  =  4-4  (2-i3) 

At  this  point,  a  brief  comment  about  notation  is  required.  In  Chapter  I,  the  symbol 
“y”  is  used  to  represent  fractional  frequency.  In  this  chapter,  however,  the  same  symbol  is 
used  along  with  sub  and  super  indices  to  represent  time  differences  between  real  clocks 
from  an  ensemble.  There  should  be  no  confusion  since  the  context  in  which  this  symbol 
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appears  clearly  identifies  its  meaning.  Moreover,  when  the  symbols  y  ,  y(t) ,  y  or  y.  are 
used,  they  always  represent  fractional  frequency  deviation.  On  the  other  hand,  when 
superscripts  are  present  as  in  yl] ,  y\ ,  y‘  or  yij  they  represent  time  differences  between 

pairs  of  real  clocks  from  the  ensemble  of  clocks.  To  avoid  confusion,  the  subscripts  R  or 
S  will  be  always  added  to  the  “y”  when  representing  fractional  frequency  deviation  as  in 
y'Rk  or  y‘Sk  respectively  when  referring  to  absolute  or  relative  fractional  frequency 
deviation. 


Now,  it  is  possible  to  think  of  yiJ  as  a  different  stochastic  process  that  relates  how 
the  time  differences  between  clocks  i  and  j  evolve  with  time.  This  process  can  be  also 
characterized  in  terms  of  its  expected  value  and  variance.  In  particular,  the  variance  of  the 
process  yj  will  be  a  function  of  the  variances  and  covariance  of  the  clocks  i  and  j  ( rn ,  r 
and  r..).  The  goal  here  is  to  estimate  the  absolute  quantities  r..  associated  with  each  clock 
by  using  relative  measurements  among  all  pairs  of  clocks. 


With  an  ensemble  of  K  clocks  it  is  possible,  in  theory,  to  obtain  X  fX  -  \)l  2 

processes,  one  per  pair  of  different  clocks.  However,  for  real  applications  only  K- 1  of 
them  are  independent.  The  reason  for  this  is  that  from  the  two  sources  of  noise  present  in 
the  model,  the  internal  noise  associated  with  each  clock  and  the  measurement  noise,  the 
last  one  is  negligible.  By  only  measuring  the  time  differences  between  the  first  K- 1 
clocks  and  the  Xth  it  is  possible  to  find  the  differences  between  each  other  pair  just  by 
taking  the  difference  of  measurements  with  respect  to  the  Xth. 


When  the  Xth  is  chosen  as  a  reference  and  the  remaining  X-\  clocks  are  compared 
to  it  at  N  different  instants  of  time,  the  set  of  time-series  measurements  can  be  arranged  in 
the  following  A-vectors:  y1A  =  xl  - xk  ,  y2k  =  x1  - xk  ...  yk~hk  =  xk~l  - xk  .  And 
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finally,  these  K- 1  processes  containing  the  measurements  can  be  arranged  in  the 
following  Nx(K- 1)  matrix: 

Y  =  (y'K,y1K,y,K,...,yK-'-K)  (2.14) 

The  expected  value  of  the  processes  y‘K  can  be  estimated  as  a  function  of  its 
elements  as: 


—iK 

y 


+  y2  +....+ 


and  arranged  into  a  constant  vector  of  N  elements  as: 


r 


(2.15) 


(2.16) 


Also,  for  purposes  of  compactness,  the  expected  values  associated  with  all  of  the 
processes  for  the  clocks  in  the  ensemble,  can  be  arranged  in  the  following  N  x  (K  - 1) 
matrix: 

rK,rK,-,r^)  (2.i7) 


If  the  K—  1  processes  of  relative  time  deviations  y,k  are  arranged  in  the  vector  v 
v  =  (y'k ,  y2k ,  yu , ...,  yk  u )  ,  then  the  covariance  matrix  associated  with  the  processes 
can  be  defined  as: 

S  =  £{vV}-£{v}-£{vr}  (2.18) 


Now,  the  covariance  matrix  of  the  processes  associated  with  the  available  samples 
for  the  K- 1  pair  of  clocks  in  the  ensemble  can  be  estimated  by  [20,21,23]: 


s=— (y-y)t\y-y) 

N- lv  ’  1  ’ 


(2.19) 


The  matrix  S  is  referred  to  as  the  relative  covariance  matrix  (in  contrast  with  the 
absolute  covariance  matrix  R)  since  it  is  the  covariance  matrix  involving  the  processes 
relating  the  time  differences  between  real  clocks  in  the  ensemble  with  respect  to  another 
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real  clock  of  the  ensemble  that  is  chosen  as  the  reference.  A  generic  element  of  the  matrix 
S  can  be  expressed  as: 

s»=jh  -y“f  ■(** -y*)  <2-20> 

Each  element  of  S  represents  a  variance  estimate  if  (/  =  j)  and  is  positive,  or  a 
covariance  estimate  if  (/  ^  j)  and  can  be  positive  or  negative. 


From  now  on,  the  superscript  K  will  be  omitted  since  the  last  clock  will  always  be 
chosen  as  the  reference. 

D.  CONSIDERATIONS  ON  THE  COVARIANCE  MATRICES  R  AND  S 


The  absolute  and  relative  covariance  matrices  R  and  S  have  been  derived  for  the 
processes  of  time  deviations  obtained  from  comparative  readings  of  clocks.  However,  the 
results  for  R  and  S  are  completely  general  for  any  other  type  of  process.  In  particular,  in 
frequency  stability  the  processes  that  are  commonly  used  are  those  obtained  with 
fractional  frequencies  instead  of  time  deviation.  Since  there  is  a  linear  relationship 
between  time  deviation  and  fractional  frequencies  (1.10)  it  is  straightforward  to  convert 
from  one  to  another  and  use  fractional  frequencies  as  the  characteristic  process.  In  this 
case,  the  matrices  R  and  S  represent  a  measure  of  variability  of  the  internal  frequency 
generated  inside  each  clock.  The  less  variability  in  the  frequency  corresponding  to  a 
particular  clock  the  more  stable  it  is  and  the  higher  its  weight  should  be. 


Our  object  is  to  find  the  absolute  covariance  matrix  R.  Unfortunately,  the 
elements  of  this  matrix  are  not  directly  observable  since  the  hypothetical  perfect  clock 
does  not  exist.  On  the  other  hand,  the  relative  covariance  matrix  of  fractional  frequencies 
S  can  be  estimated  from  relative  measures  between  clocks  from  the  ensemble.  The  next 
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section  focuses  on  the  relationship  between  these  two  matrices  and  the  mathematical 
tools  needed  to  obtain  information  about  R  starting  from  the  estimation  of  S. 

1.  Relationship  Between  the  Matrices  R  and  S. 


The  absolute  and  relative  covariance  matrices  R  and  S  are  very  closely  related  and 
can  be  estimated  as  functions  of  the  time  comparison  of  the  K- 1  pairs  of  clocks  that  form 
the  ensemble.  The  dimension  of  the  matrix  R  is  KxK  while  the  corresponding 
dimension  for  S  is  (A'-l)x(A'-l).  Both  are  symmetric  and  positive  scmidefinitc, 

consequently,  there  are  K-(K  + 1 ) / 2  independent  elements  in  R  and  (K  -  K /2 

independent  elements  in  S.  The  matrix  S  can  be  estimated  directly  from  the  time 
measurements  of  the  pairs  of  clocks,  but  the  matrix  R  has  K  more  elements  than  S  so 
there  is  not  a  unique  determination  of  R  and  an  optimization  criteria  should  be  defined  in 
order  to  find  the  “best”  matrix  R. 


The  relationship  between  R  and  S  can  be  derived  by  observing  that  the 
measurement  matrices  X  and  Y  are  related  by: 

Y  =  X  H,  (2.21) 

where  H  is  the  K  x  (K  - 1)  matrix: 


H  = 


(2.22) 


being  /  the  (K  - 1)  x  (K  - 1)  identity  matrix  and  uT  the  ( K  - 1)  vector  of  ones: 

«  =  (1,1,1,...,  if  (2.23) 

The  matrices  containing  the  expected  values,  X  and  Y ,  are  also  equally  related 
with  matrix  H  as:  Y  =  X  ■  H 
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So  that,  finally,  we  have  the  following  fundamental  relationship  between  matrices 
R  and  S: 


S  =  H  R-  H 


(2.24) 


For  reasons  that  will  become  clear  later,  the  covariance  matrix  R  is  partitioned 
separating  the  free  and  non-free  parameters  according  to  the  following  structure: 


R  = 


R 


\FK- 1 


(2.25) 


'KK  ) 


Here  the  matrix  R  is  of  order  ( /f  -  1 )  x  ( /f  -  1 )  and  contains  the  first  K- 1  rows  and 

columns  of  R,  and  vector  rK]  is  given  by  rK_x  =  {rlK,  r2K,  riK,...,  rK_lK^  .  The  relationship 
between  R  and  S  can  be  stated  as: 


R 


S  =  HT  R  H  =  [I  -u) 

R  +  rKK-\_u  uT~^-u-rK^T  -r^  u 


\rK- 1 


KK  J  L 


-u 


(2.26) 


This  last  equation  shows  that  R  can  be  uniquely  obtained  from  the  matrix  S  once 
the  elements  contained  in  vector  rK1  and  scalar  rKK  are  known. 


During  the  process  of  estimation  of  clock  instabilities,  measurement  data  are 
usually  pre-filtered  to  use  for  example,  the  Allan  variance  statistic  instead  of  the  classical 
variance/covariance.  The  derivation  of  matrices  R  and  S  is  completely  general  and  the 
relationship  between  them  also  holds  when  Allan  variances  and  covariances  are  used 
inside  S  [21,  25].  This  is  the  approach  followed  in  this  thesis  and,  unless  otherwise  stated, 
the  variances  and  covariances  used  here  will  be  the  Allan  (co)variances  associated  with 
the  clock  ensemble. 

In  order  to  obtain  the  matrix  S,  the  Kth  clock  of  the  ensemble  is  chosen  as  a 
reference  and  all  of  the  pairwise  comparisons  among  the  other  K  —  1  clocks  and  the  K,h 
are  used  in  the  computation  of  S. 
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2.  A  Simple  Case  with  an  Ensemble  of  K= 3  clocks. 

An  example  of  an  ensemble  of  three  clocks  is  presented  here  for  several  reasons: 

•  It  is  possible  to  derive  explicit  expressions  for  the  elements  of  R  as 
functions  of  the  corresponding  elements  of  S  and  the  problem  is 
mathematically  tractable  [20]. 

•  In  the  more  general  case  when  K  >  3  the  same  approach  will  be 
followed  but  a  final  solution  for  matrix  R  will  be  obtained  by 
numerically  solving  an  optimization  problem  [21,22], 

•  There  is  an  easy  transition  to  the  classical  and  widely  used  Three- 
Cornered  Hat  method  [23], 


If  the  number  of  clocks  is  three  and  the  third  is  chosen  as  a  reference,  the 
available  data  are  the  N  comparison  samples  between  the  pair  of  clocks  1  and  3  and 
between  2  and  3.  The  elements  of  the  matrices  R,  S  and  H  are  then: 


(2.27) 


(2.28) 


'  1  0^ 

H=  0  1  (2.29) 

-1  ~b 

And  by  using  the  relationship  between  them  (2.24): 

■hi  =  h i  -  2  •  rn  +  r33  s12  =  rn  -  ru  -  r23  +  r33 ^ 

S12  =  r\2  ~~  ri3  ~  r23  r33  ^22  =  ^22  _  ^  '  ^23  +  ?33  J 


S  =  Ht  R  H  = 
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There  are  3  equations  in  (2.30)  corresponding  to  su,  s22)  that  involve  6 
unknown  parameters:  (rM,  rl2,  rl3,  r22,  r23,  r33  )  ,  consequently,  there  are  infinitely  many 
solutions  to  the  problem  of  finding  them.  If  the  elements  of  R  associated  with  the 
reference  clock  (r13,  r23,  r33)  are  considered  as  free  variables,  then  the  remaining  elements 

of  R  can  be  expressed  as  functions  of  them  and  the  elements  of  S  and  the  three  equations 
become: 

rn=sn-r33+2-ri3 

rn=sn-r33+r]3+r23  (2.31) 

r22  =  S22  ~  ^33  2  '  f>3 

However,  not  all  of  the  solutions  are  feasible.  There  are  constraints  to  satisfy 
since  the  matrix  R  must  be  positive  semidefinite.  For  a  general  matrix  to  be  positive 
semidefinite  it  is  necessary  and  sufficient  that  all  the  detenninants  of  the  principal  minors 
of  the  matrix  are  positive.  It  could  be  thought  that  there  are  as  many  constraints  as  the 
number  of  principal  minors  in  the  matrix  R ,  that  is,  3,  but  it  is  not  the  case  here  because 
the  matrix  R  is  not  a  general  matrix.  In  [20,21]  it  is  shown  that  since  matrix  R  is  obtained 
from  matrix  S  which  is  already  positive  semidefinite,  the  only  necessary  and  sufficient 
conditions  for  R  to  be  positive  semidefinite  is  that  the  detenninant  of  R  must  be  non¬ 
negative: 

\R\  =  >\  I’  r22  ■  r32  +  2  '  ri2  ■  >23  '  f  3  ~  4  '  r22  ~  f  2  '  ^33  “  4  ‘  >\  1  ^  0  (2-32) 

This  means  that  there  are  now  3  equations  (2.31),  1  constraint  (2.32)  and  6  unknowns  . 
There  are  still  infinitely  many  solutions. 

3.  Selection  of  the  Free  (Co)variances  r13,  r23  and  r33 . 


An  optimality  criterion  must  be  defined  to  reach  a  solution  for  the  elements  of  the 
covariance  matrix  R.  The  objective  function  proposed  by  [20]  includes  a  global  measure 

of  covariance  of  the  ensemble:  quadratic  mean  covariance  G(R)  =  ^(r13  +r13  +r23)/3  , 
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which  can  be  expressed  as  a  function  of  the  free  parameters  grouped  as 
r3  =  (r13 ,  r23 ,  r23  )7  ,  and  the  elements  of  the  relative  covariance  matrix  S  by  using  (2.30): 


[G(R)f  ;+4)/3  = 


[G(r3;5)]2  = 

(  2  •  (f3  -  ^33  )2  +  2  •  {rn  -  r33 )  •  (fi3  -  )  +  -  (2.33) 

2  ■  [r23  —  /33  )  +  2  •  (2  •  r33  +  sl2 )  •  (rn  —  r33 )  + ... 

2  '  (2  ‘  r33  +  S\2  )  '  (r23  _  r33  )  "*"  2  '  r33  +  (512  +  ^33  ) 


Equation  (2.33)  can  be  generalized  to  any  number  of  clocks  by  applying  the  basic 
linear  relationship  S  =  HT  ■  R-  H  in  order  to  express  G(R )  as  a  function  of  the  free 

parameters  rK  and  the  elements  of  S:  \G(rK;  51)]" .  It  will  be  used  as  a  component  of  two 

different  objective  functions  that  will  be  analyzed  below.  A.  Premoli  and  P.  Tavella  [20] 
suggest  an  unconstrained  optimization  problem  to  find  r13 ,  r23  and  r33 .  They  minimize  by 

choice  of  rl3,  r23  and  r33  the  following  objective  function: 


F(r3'S)  = 


\R\ 


(2.34) 


The  idea  is  to  minimize  the  “global  correlation”  among  clocks  while  maintaining 
the  positive  semidefiniteness  of  R.  The  combination  of  G(r3;  S )  in  the  numerator  and  \R\ 
in  the  denominator  finds  the  matrix  R  with  smallest  global  correlation  and  at  the  same 
time  ensures  the  positive  semidefiniteness  of  R  by  the  presence  of  \R\ .  If  during  the 
minimization  process  the  incumbent  solution  approaches  the  boundary  of  the  feasible 
region  |i?|  =  0,  then  the  objective  function’s  value  increases  showing  that  this  is  not  a 
valid  direction  to  minimize. 


For  the  case  of  3  clocks  A.  Premoli  and  P.  Tavella  showed  that  the  solution  to  the 
minimization  problem  is  unique  and  can  be  performed  analytically  to  find  the  optimal 
values  r*3 ,  r23  and  r33 . 
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4. 


The  Classical  Three-Cornered  Hat  Method. 


The  classical  method  most  widely  used  to  find  the  variances  associated  with  each 
clock:  rn,  r22  and  r33  can  be  derived  as  a  particular  case  of  the  method  proposed  above. 
The  three-cornered  hat  method  [20,26]  assumes  that  the  clocks  are  completely 
independent  and  uncorrelated.  It  means  that,  for  the  same  case  of  three  clocks, 
rn  =  r13  =  r23  =  0 .  Under  this  hypothesis,  since  the  covariance  matrices  R  and  S  are  related 

by  S  =  Ht  R  H ,  equation  (2.30)  becomes: 


5  = 


^11  =hl+r33 


r33 


s22  r22  +  r33  j 


(2.35) 


In  order  for  R  to  be  positive  semi  definite,  since  rn,r22,r33  >0  and  also  su,s22>  0,  the 
following  conditions  must  be  met: 

s12>0,  sl2Zsn,  sn<s22  (2.36) 


In  that  case,  the  solutions  for  rn,  r22  and  r33  are  the  same  that  the  ones  obtained  in  the 
minimization  problem  and  have  a  simple  representation  obtained  from  (2.3 1): 

rn  =sn-sn 

r22  —  S22  ~  ^12  (2.37) 

r22  =  S\2 

Since  the  solution  involves  a  system  of  three  simultaneous  equations  with  three 
unknowns;  there  is  no  need  to  perfonn  a  minimization  since  there  are  no  free  parameters  , 
and  the  condition  of  positive  definiteness  is  consequently  not  imposed.  As  a  consequence, 
the  solution  sometimes  may  produce  negative  values  for  the  variances,  which  is  not 
acceptable. 


Several  approaches  to  avoid  the  appearance  of  negative  variances  have  been 
proposed,  such  as  considering  only  the  absolute  value  [27]  or  substituting  its  value  by 
zero  when  negative  results  are  obtained.  None  of  them  is  based  on  theoretical  or  physical 


37 


considerations.  They  are  just  a  means  to  avoid  this  nuisance.  The  causes  leading  to  the 
estimation  of  negative  variances  when  the  popular  three-cornered  hat  method  is  used  are: 

•  Uncertainty  in  the  measured  time  differences. 

•  Lack  of  contemporaneity  of  measurements. 

•  Presence  of  correlation  between  clocks. 


When  some  of  the  above  conditions  are  not  met,  the  three-cornered  hat  method  is 
no  longer  applicable.  The  introduction  of  free  parameters  requires  the  definition  of  a 
suitable  minimization  criterion  in  order  to  provide  the  analytical  conditions  to  obtain  the 
desired  values  for  the  free  parameters.  The  optimality  criterion  is  usually  based  on  the 
minimization  of  some  function  that  measures  a  global  covariance  associated  with  the 
whole  ensemble  of  clocks. 

5.  Optimal  Selection  of  the  Free  Covariances:  Constrained  or 
Unconstrained  Minimization? 


When  the  number  of  clocks  in  the  ensemble  is  greater  than  three,  it  is  shown  in 
section  C.l  that  the  problem  of  finding  the  covariance  matrix  R  from  the  covariance 
matrix  of  observations  S  is  under  determined.  In  particular,  if  the  number  of  clocks  in  the 
ensemble  is  K,  there  are  K-(K  + 1 )  /  2  independent  elements  in  R  and  if-(W-l)/2 

independent  elements  in  S;  therefore  there  are  always  exactly  K  free  parameters  to 
detennine  by  means  of  a  minimization  problem.  Unless  otherwise  stated,  in  the  rest  of 
this  paper  the  election  of  the  free  parameters  will  be  (co)variances  associated  with  the 
clock  that  is  chosen  as  the  reference  (which  is  also  usually  chosen  as  the  one  that  is 
believed  to  have  the  best  stability  characteristics),  i.e.  the  Xth  clock  in  the  ensemble.  Then 
r\K’r2K’r3K’—’rKK  constitute  the  decision  variables  of  the  problem.  For  compactness, 

group  these  in  the  vector  rK  =  (r]K,  r2K,  nK,  ...,rKK)T  .  The  rest  of  the  elements  of  the  matrix 
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R,  that  is,  the  ones  that  do  not  involve  the  Xth  clock  can  be  always  expressed  as  functions 
of  the  elements  of  S  and  the  decision  variables: 


=  fy(rK’S)  Vi  Je  {1,2, 3,...,*-!} 


(2.38) 


The  approaches  followed  in  [21]  and  [25]  suggest  solving  an  unconstrained 
minimization  problem  with  an  objective  function  taking  one  of  the  following  fonns: 


where 


Ft(rt;S) 

f2{ks) 


Kv'bl 

|*K;s)f 


(2.39) 


[c(vS)MZ''»! 


j= 1  ><j 


(2.40) 


The  presence  of  the  square  of  the  determinant  of  R  in  the  denominator  provides, 
according  to  Torcasso  [25],  a  better  numerical  convergence  and  helps  to  maintain  the 

incumbent  solution  inside  the  feasible  region  given  by  | R  (  0,  In  [21]  it  is  shown 

that  this  expression  has  a  nice  closed  form  expression  suitable  for  numerical 
computations: 


|*^  )  '  yKK  ~  \rK- 1  _  rKK  '  u\  j  X  S  X  (r KK  ~  ,  —  rKK  •  lljj 


(2.41) 


6.  New  Proposal  for  the  Constrained  Minimization  Problem. 


The  solution  adopted  in  this  paper  is  based  on  a  generalization  of  the  solution 
proposed  in  [22,24]  which  has  an  immediate  physical  interpretation.  As  a  side  effect,  it 
can  be  used  on  clocks  in  an  ensemble  consisting  of  clocks  from  different  laboratories. 
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The  objective  function  based  on  the  quadratic  mean  covariance 

2  K  K 

[G(rA.;5')]  has  the  attractive  property  of  convexity  in  the  space  of  the  free 

i=l  i< j 

parameters  rK.  However,  the  (co)variance  tenns  are  not  normalized;  they  can  take  any 
value  and  can  oscillate  orders  of  magnitude  from  clock  to  clock  -  especially  if  the  clocks 
do  not  have  the  same  characteristics.  The  final  solution  can  be  greatly  influenced  by  the 
characteristics  of  one  particular  clock.  If  the  objective  function  is  modified  to  be  a 
“quadratic  mean  correlation”  of  the  ensemble,  instead  of  covariance,  then  the  effect  of 
one  clock  presenting  big  variance  and  covariance  terms  does  not  have  such  a  big  impact 
in  the  final  solution  obtained  by  the  minimization.  This  effect  will  be  shown  with  an 
example.  This  new  objective  function,  the  one  that  will  be  minimized  in  this  thesis,  has 
the  following  expression: 


F,(> 


r 1 

y  y  '  ij 


7=1  i<j 


■rjj 


(2.42) 


In  the  approach  presented  here,  one  new  set  of  constraints  has  been  added  to  the 
minimization  problem.  To  the  knowledge  of  the  author,  no  reports  accounting  for 
negative  correlations  between  pairs  of  clocks  under  similar  environmental  conditions 
have  been  reported,  so  it  makes  sense  to  constrain  the  covariances  or  correlation 
coefficients  to  only  non-negative  values.  Therefore,  clocks  from  the  same  laboratory  will 
be  constrained  to  present  positive  correlations  between  them.  Correlation  between 
laboratories  is  unrestricted  and  consequently  can  take  both  positive  or  negative  values. 
Negative  correlations  between  laboratories  could  be  caused  by  diurnal,  seasonal  or 
geographical  factors. 


7,  Mathematical  Formulation  of  the  Minimization  Problem. 


The  general  minimization  problem  to  be  solved  when  two  different  laboratories 
are  involved  is: 
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(2.43) 


Min 

rK 

s.t. 


j= i  «</ 


|/?(r,;5)|>0 

r7(rA,;5)>0  V /,  y  e  Cj,  /  <  y 
r./(fAy4S)>0  V  i,  j  eC2,i  <  j 


Where  C,  contains  the  indices  of  all  clocks  belonging  to  laboratory  1  and  C2 
contains  the  indices  of  all  clocks  belonging  to  laboratory  2.  The  requirement  of  positive 
correlation  among  the  clocks  in  the  ensemble  is  assured  by  the  constraints  rJrK;S)  >  0  . 


For  the  simple  case  of  an  ensemble  of  K  clocks,  all  from  the  same  laboratory,  the 

K-IK  + 1) 

absolute  covariance  matrix  „  consists  of  - - - -  different  elements 


ru ,  i 


i,j  e  {1,2,..., if}  that  need  to  be  detennined.  The  relative  covariance  matrix 


K-(K-l) 

S,„  ,w„  has  - - - -  different  elements  that  are  known  from  measured  data.  The 

(x-i)x(x-i)  2 

K  •( K  —  W 

relationship  (2.24)  among  R  and  H  then  provides  - - -  distinct  linear  equations 


that  relate  the  elements  r..  and  su ,  so  there  are  only 


^•(^  +  1)  K-(K-l) 


=  K  free 


parameters  among  the  elements  of  R.  They  are  the  decision  variables  and  will  be  chosen 
as  the  elements  of  R  associated  with  the  Xth  clock:  rK .  Before  carrying  out  the 
minimization  problem  it  is  necessary  to  express  the  objective  function  and  the  constraints 
as  functions  of  only  the  decision  variables  rK .  Then,  once  the  minimization  problem  has 

been  solved,  the  remaining  elements  of  R  can  be  uniquely  computed  using  again  the 
linear  equations  that  relate  the  matrices  R ,  H  and  S. 
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E. 


ESTIMATION  OF  THE  COVARIANCE  MATRICES  R  AND  S 


1.  Notation  and  Definitions. 

Let  the  quantity  x'm  be  a  realization  of  the  process  x'  that  represents  the  time 
difference  between  clock  i  and  an  ideal  clock  at  the  instant  of  time  in  To,  where  To  is  the 
sampling  interval.  Similarly,  let  y'm  be  a  realization  of  the  process  y  that  represents  the 
time  difference  between  clock  i  and  the  Xth  clock  or  reference  at  the  instance  of  time  m  To. 

According  to  these  conventions,  the  absolute  and  relative  fractional  frequencies  can  be 
expressed  using  (1.11)  and  the  new  notation  as: 


yiRk=yiM  =  h±^L 

(2.44) 

And 

ro 

(2.45) 

For  k  >  1 ,  for  intervals  of  time  multiples  of  r0 ,  the  usual  average 

of  the  fractional 

frequency  (1.19)  is  always  used: 

Y  m 

fSk(T  =  m  ■  t0)  =  —  )»+/ 

(2.46) 

m  /=1 

Y  m 

y'Rk(T  =  m  ■t0)  =  ~Yj  ylR{k~\)m+i 

m  i= i 

(2.47) 

Recall  that  both  x‘  and  y‘  represent  the  differences  between  clock  i  and  an  ideal 
clock  and  between  clock  i  and  clock  K.  Consequently,  y‘Rk  and  y‘Sk  represent  average 
variations.  If  clock  i  is  perfect  then,  y'Rk  will  be  zero,  but  y‘Sk  does  not  have  to  be  zero 
since  clock  K  may  not  be  perfect.  If  clock  i  is  not  perfect  then  y‘Kk  will  characterize  clock 
i,  but  y'Sk  will  include  the  instabilities  associated  with  both  clocks  i  and  K.  For  each  of 
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the  first  K- 1  clocks  in  the  ensemble,  y'sl(  represents  a  time  series  that  will  be  used  to 
estimate  y'Rk . 

2.  Estimation  of  Allan  (Co)variances  from  the  Measured  Data. 


The  fractional  frequency  deviation  time  series  associated  with  the  K- 1  pairs  of 
clocks  could  be  directly  used  to  compute  the  (A-l)x(W-l)  covariance  matrix  S  using 

the  raw  data.  However,  this  approach  is  never  used  in  practice  since  there  may  be  present 
different  types  of  noise  for  which  the  usual  variance  does  not  converge  as  previously 
shown  in  Figure  2.  Instead,  the  raw  data  are  filtered  in  the  time  domain  to  obtain  the 
covariance  matrix  S  based  on  the  Allan  covariances  instead  of  the  traditional 
(co)variances. 


Using  the  previous  notation,  the  Allan  variance  associated  with  clock  i  referenced 
to  clock  K  can  be  expressed  using  ( 1 . 13)  as: 

5,7  (r>  =  |((^+i(r>_^(r))2}  (2-48) 

and  the  Allan  variance  associated  with  clock  i  referenced  to  the  ideal  clock  can  be 
expressed  as: 

r«(T)  =  ^((kLiW-KfA))2)  (2.49) 

where  again,  the  (  denotes  mathematical  expectation  or  infinite  time  average.  We 

assume  that  the  process  of  averaged  fractional  frequencies  is  stationary  and  ergodic  so 
that  the  last  equation  is  well-defined. 
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Similarly,  when  (/'  ^  j)  the  Allan  covariances  of  clocks  i  and  j  referenced  to 

clock  K  or  referred  to  the  ideal  clock  also  can  be  written  respectively,  according  to  (1 . 15) 
as: 

Sij(r)  =  ^((yLAT)-yL(T)){yLi(T)-yL(T)))  (2-5°) 

rij(T)  =  |  {(ylRkJT)  -  yRk  0))  •  (yLM  -  yi*  (*■)))  (2-5 1) 

As  usual,  the  covariance  reduces  to  variance  when  (/'  =  j). 


The  covariance  matrix  S  can  now  be  directly  estimated  from  the  filtered  data 
obtained  from  the  time  comparisons  of  the  K  clocks.  The  overlapped  estimator  given  by 
[3]  for  a  given  r  =  m  •  r0  is: 


s'(r=,"'r">==r^) 

N-2n 


(2.52) 


X  {y'sk+\ (m  ■  h)  - y'sk (fn  •  ro))  •  (yLi (m  •  h)  - yk (m  ■  h)) 


k= 1 


and  similarly  for  the  elements  of  the  covariance  matrix  R: 

1 


rij(j  =  m  ■  r0)  = 

N-2n 


2{N-2m) 


(2.53) 


X  (k/«+1  (m  ■  T0)  - y'Rk  (m  •  r0))  •  (yJRk+l (m  ■  r0)  -  yJRk  (m  ■  r0)) 


k= 1 


The  elements  of  the  matrix  S  are  computed  from  the  time  deviation  series  using 
equations  (2.52),  (2.46)  and  (2.45). 
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III.  COMPARISON,  TESTING  AND  SIMULATION 


A.  OVERVIEW 

1.  Description  of  the  Different  Tests  Performed 


The  theoretical  formulation  of  the  method  to  compute  the  variability  of  the 
fractional  frequency  deviation  associated  with  each  particular  clock  in  the  ensemble  has 
been  stated  in  the  last  chapter.  The  next  step  in  the  analysis  of  the  validity  of  the  method 
involves  some  comparisons  with  results  taken  from  other  researchers,  and  also  with 
simulations  to  test  the  accuracy  and  the  response  to  dynamic  changes  in  the  parameters  of 
the  clocks. 


Four  different  tests,  increasing  in  complexity,  are  summarized  below: 

•  A  point  comparison  of  the  method  for  one  particular  artificial  case 
analyzed  by  several  authors.  Comparable  results  with  what  others  have 
obtained  will  give  some  insight  about  the  possibilities  of  the  new  method. 

•  A  simulation  involving  the  generation  of  random  covariance  matrices  as 
the  starting  point  and  use  the  5  different  methods  described  below  in 
section  III  B,  to  solve  the  same  problem.  Three  different  MOP' s  will  be 
defined  and  used  to  compare  the  results  provided  by  each  approach. 

•  A  simulation  to  test  the  dynamic  behavior  of  the  proposed  method.  Instead 
of  generating  random  absolute  covariance  matrices  R ,  the  variance  tenns 
associated  with  each  clock  will  be  changed  according  to  a  predefined 
function  to  test  the  sensitivity  against  small  changes  in  the  input  data. 

•  The  last  simulation  involves  the  solution  of  a  complete  problem.  The 
phase  or  time  data  of  each  clock  is  generated  including  additive  noise  with 
predefined  spectral  characteristics  and  Allan  variances.  This  simulation 
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involves  the  computation  of  the  relative  covariance  matrix  S  from  phase 
data,  the  minimization  problem  to  derive  the  absolute  covariance  matrix  R 
and  the  weights  of  each  clock. 

2.  Analysis  of  the  Uniform  Random  Number  Generator. 


Before  entering  into  the  specific  details  of  the  simulation  it  is  worth  saying  that 
the  software  package  used  is  matlab©  (version  6.5)  with  the  optimization  toolbox 
(version  2.2).  Extensive  use  of  the  Uniform  Random  Number  Generator  (URNG)  built 
into  matlab  is  required.  Unfortunately,  there  is  no  information  available  about  this 
particular  URNG  and  so,  before  using  it,  the  tests  included  in  the  diehard  suite  of  tests 
created  by  Marsaglia  [28]  are  applied  to  a  24MB  file  containing  uniform  random 
numbers  generated  in  matlab.  The  steps  to  test  the  generator  are  summarized  below: 

•  Generate  32-bit  unsigned  integers  by  mapping  the  range  of  the  interval  of 
the  URNG  [0, 1]  into  the  integer  interval  [0, 2 '2  - 1]  as  required  by  the  test 
specifications  [28], 

•  Conversion  of  the  32-bit  integers  into  8-character  hexadecimal  numbers. 

•  Create  a  24MB  text  file  consisting  of  80  character  lines  (10  numbers)  with 
the  numbers  previously  generated. 

•  Create  a  binary  file  using  the  last  data  file  as  input 

•  Finally,  execute  the  test  program  and  provide  the  binary  file  as  input. 

The  results  indicated  that  the  URNG  is  of  sufficient  quality  to  proceed 

3.  Definition  of  Proposed  Measures  of  Performance  (  MOP  ). 


In  order  to  test  and  compare  the  results  derived  in  the  simulation,  three  measures 
of  performance  are  defined.  They  are  based  on  the  absolute  covariance  matrix,  R ,  and  its 


46 


correlation  matrix,  p .  Just  by  looking  at  Table  9,  is  difficult  to  say  which  of  them 
provide  results  closer  to  the  true  values  of  R  and  p .  In  order  to  compare  them,  three 
different  measures  of  performance  {MOP)  are  defined: 

MOP] :  Sum  of  the  squared  residuals  of  the  elements  of  the  computed  R. 

M°p>=YL{r«-fif  <31> 

7=1  tej 

MOP2 :  Sum  of  the  squared  residuals  of  the  elements  of  the  correlation  matrix  p  . 

MOP2=t,t,(p.,-Plf  O'2) 

7=1  tej 

MOP, :  Sum  of  the  squared  residuals  of  the  elements  corresponding  to  the 
variances  of  R,  that  is,  the  diagonal  terms  of  R. 

MOP,  =  f.(, f  (3.3) 

M 

The  third  MOP  is  the  simplest  and  the  preferred  one;  the  most  important  values  in 
the  matrix  R  are  the  variance  terms  located  in  the  main  diagonal,  since  they  are  directly 
used  to  compute  the  weighting  factors  for  each  clock. 


B.  TEST  N°l,  SINGLE  POINT  COMPARISON 


In  order  to  compare  the  approach  presented  here  with  previous  approaches,  the 
same  sample  matrix  S  that  has  been  used  by  other  authors  [22,24],  is  used  below  so  that 
the  final  solutions  can  be  easily  compared.  The  five  methods  used  in  the  analysis  are 
listed  below: 

•  1 .  Three-cornered  hat  method  [20] 

•  2.  Tavella  method  [20] 
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•  3.  Torcasso  et  al.  method  [25] 

•  4.  Galindo  F.  J.  et  al.  method  [22,24] 

•  5.  Proposed  Ruiz-Galindo  method. 


The  starting  (simulated)  covariance  matrix  R  is  based  on  an  ensemble  of  four 
clocks  having  different  characteristics  and  also  includes  correlation  between  different 
pairs  of  clocks.  The  covariance  matrix  R  for  this  example  is  given  below.  According  to 
[22],  it  is  generated  with  a  correlation  coefficient  for  each  pairs  of  clocks  fixed  at 

Pij  =0.10 : 


>i 

rn 

6  3 

r  > 
'14 

"2.085 

0.153 

2.237 

0.297 

>h 

r22 

r2A 

0.153 

1.143 

1.744 

0.234 

rn 

2.237 

1.744 

239.62 

3.524 

Jia 

r2A 

v0.297 

0.234 

3.524 

4.064 

The  computed  correlation  coefficients  associated  with  R  are  close  to  0.10: 


A  2 

As 

a4^ 

"1.000 

0.099 

0.100 

0.102" 

A  2 

P22 

Pi. 3 

PlA 

0.099 

1.000 

0.105 

0.109 

P\2 

Pi. 3 

Pi3 

Pm 

0.100 

0.105 

1.000 

0.113 

\P\A 

PlA 

Pm 

PaaJ 

v0.102 

0.109 

0.113 

1.000, 

The  covariance  matrix  S  corresponding  to  R  can  be  computed  using  the 
relationship  between  R  and  S  : 


(  sn  si2 

S13  'l 

f  5.555 

3.686 

2.480  " 

S  =  Hr  R  H  = 

S12 

S22 

S23 

= 

3.686 

4.739 

2.050 

S23 

S33  ) 

^2.480 

2.050 

236.636, 

where 

'+1  0  0" 
0+10 

H  = 

0  0+1 

v-1  “I  "I, 
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Now,  using  the  5  different  methods,  the  results  for  R  and  p  are  summarized  in 
Tables  9  and  10  below: 


Method 

R 

P 

True 

Value 

R,p 

'  2.085  0.153  2.237  0.297^ 

0.153  1.143  1.744  0.234 

2.237  1.744  239.62  3.524 

v  0.297  0.234  3.524  4.064 , 

'  1.000  0.099  0.100  0.102^ 

0.099  1.000  0.105  0.109 

0.100  0.105  1.000  0.113 

v  0.102  0.109  0.113  1.000, 

Method  1 

Three 

Cornered 

Hat 

r  2.817  0.000  0.000  0.000 " 

0.000  2.000  0.000  0.000 

0.000  0.000  233.90  0.000 

v  0.000  0.000  0.000  2.739, 

'  1.000  0.000  0.000  0.000" 

0.000  1.000  0.000  0.000 

0.000  0.000  1.000  0.000 

v  0.000  0.000  0.000  1.000, 

Method  2 

(Tavella) 

r  2.265  0.625  -0.134  -0.205" 

0.625  1.907  -0.334  0.024 

-0.134  -0.334  234.70  0.471 

v -0.205  0.024  0.471  2.880 , 

'  1.000  0.301  0.006  0.080^ 

0.301  1.000  0.016  0.010 

0.006  0.016  1.000  0.018 

v  0.080  0.010  0.018  1.000, 

Method  3 

(Torcasso) 

r  2.435  0.805  -0.135  -0.035" 

0.805  2.098  -0.326  0.205 

-0.135  -0.326  234.53  0.470 

v -0.035  0.205  0.470  3.051 , 

'  1.000  0.356  -0.006  -0.013" 

0.356  1.000  -0.015  0.081 

-0.006  -0.015  1.000  0.018 

v -0.013  0.081  0.018  1.000, 

Method  4 

(Galindo) 

'  1.859  0.003  0.035  -0.01 0 

0.003  1.069  -0.382  0.003 

0.035  -0.382  235.44  1.241 

v -0.011  0.003  1.241  3.675 ) 

"  1.000  0.002  0.002  -0.004" 

0.002  1.000  -0.024  0.001 

0.002  -0.024  1.000  0.042 

v -0.004  0.001  0.042  1.000, 

Method  5 

(Ruiz) 

'  1.867  0.003  0.425  0.000" 

0.003  1.062  0.000  0.006 

0.425  0.000  236.21  1.633 

v  0.000  0.006  1.633  3.688, 

'  1.000  0.002  0.020  0.000^ 

0.002  1.000  0.000  0.003 

0.020  0.000  1.000  0.055 

v  0.000  0.003  0.055  1.000, 

Table  9.  Comparison  of  the  covariance  and  correlation  matrices  for  the  five  approaches 
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Just  by  looking  at  the  table,  is  difficult  to  say  which  approach  provides  results 
closer  to  the  true  values  of  R  and  p  .  In  order  to  compare  them,  the  proposed  MOP's  are 
presented  in  Table  10: 


MOPx 

MOP2 

MOP » 

Method  1 

56.4062 

0.0659 

35.7756 

Method  2 

46.0122 

0.1185 

26.2352 

Method  3 

47.7893 

0 .1148 

28.0140 

Method  4 

32.4091 

0.0636 

17 . 6589 

Method  5 

21.8667 

0.0518 

11 . 8007 

Table  10.  Comparison  of  the  MOP' s  obtained  by  each  method 


For  all  MOP's,  a  smaller  value  indicates  a  better  result.  The  fifth  method  (the  new 
one)  seems  to  provide  better  results  for  all  MOP' s,  i.e.,  it  dominates  regardless  of  MOP. 

C.  TEST  N°2,  RANDOM  COVARIANCE  MATRIX  SIMULATION 

In  order  to  check  if  the  results  in  the  previous  section  are  just  a  lucky  outcome 
coming  from  a  particular  favorable  set  of  data,  a  simulation  of  similar  problems  is  made 
in  this  section. 


In  order  to  reduce  the  variance  of  the  results,  common  random  numbers  are  used 
among  the  5  different  methods.  The  simulation  is  performed  for  an  ensemble  of  4 
correlated  clocks  presenting  similar  characteristics.  It  consists  of  the  following  phases: 
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•  Generation  of  the  absolute  covariance  matrices  if4x4 . 

•  Computation  of  the  relative  covariance  matrices  S3x3 . 

•  Estimation  of  the  absolute  covariance  matrices  i?, ,  1?2 ,R3,  i?4  and  R5  from 
S  using  the  5  different  approaches. 

•  Computation  of  the  measures  of  performance  MOPl ,  MOP2 ,  and  MOP3 

•  Analyze  and  present  the  results. 


The  random  parameters  used  in  the  simulation  correspond  to  the  ones  expected  in 
an  ensemble  of  atomic  clocks,  but  with  special  attention  to  the  correlation  among  them 
since  this  is  one  of  the  objectives  of  the  present  work.  The  absolute  covariance  matrix  R 
is  obtained  by  generating  the  variance  terms  (main  diagonal)  and  desired  correlation 
coefficients  between  each  pair  of  clocks.  A  total  of  40000  runs  of  simulations  are 
performed  with  the  following  distribution  for  the  random  parameters: 

•  rn,  r22,  r33  and  rAA  uniformly  distributed  between  1  and  20  in  arbitrary 
units. 

•  pn  ,  pl3 ,  p14 ,  p23 ,  p2A  and  p3A  uniformly  distributed  between  0  and  0.1 


The  results  for  the  3  MOP ’s  are  shown  in  Figures  4,  5  and  6.  They  represent  the 
histograms  corresponding  to  the  distribution  of  each  MOP  for  the  different  approaches 
tested. 
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Histogram  MOP 1  vs.  method 


MOP 1  range 


Figure  4.  Comparison  of  1st  MOP 


Histogram  MOP2  vs.  method 
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5000 


0  0.005  0.01  0.015  0.02  0.025  0.03  0.035  0.04  0.045  0.05 


5000 


co 

4fc 


0  0.005  0.01  0.015  0.02  0.025  0.03  0.035  0.04  0.045  0.05 


5000 


3 


0  0.005  0.01  0.015  0.02  0.025  0.03  0.035  0.04  0.045  0.05 


5000 


Ln 

4fc 


0  0.005  0.01  0.015  0.02  0.025  0.03  0.035  0.04  0.045  0.05 

MOP2  range 


Figure  5.  Comparison  of  2nd  MOP 
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Histogram  MOP^  vs.  method 


MOP3  range 


Figure  6.  Comparison  of  3ld  MOP 


The  mean  and  standard  deviation  obtained  for  the  tree  MOP’s  is  shown  in  Tables 
11,  12  and  13: 


MOP  1 

Method 

#1 

#2 

#3 

#4 

#5 

Mean,  X 

4.4080 

3.8508 

3.8239 

3.5891 

2 . 0172 

Std.  Dev,  s 

3.0496 

2.7553 

2 .7427 

2 . 6750 

1.9999 

Table 

1 .  Summary  table  for  MOP  1 
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MOP's  2 

Method 

#1 

#2 

#3 

#4 

#5 

Mean,  X 

0.0200 

0.0199 

0.0182 

0.0181 

0.0104 

Std.  Dev,  s 

0.0073 

0.0323 

0.0080 

0.0072 

0.0065 

Table  i 

2.  Summary  table  for  MOP  2 

MOP  3 

Method 

#1 

#2 

#3 

#4 

#5 

Mean,  X 

2.2315 

1 . 8788 

1 .8819 

1 .7011 

1.0165 

Std.  Dev,  s 

1 . 6399 

1.3465 

1.3520 

1.2708 

0.9856 

Table  i 

3 .  Summary  table  for  MOP  3 

The  three  tables  show  that  the  proposed  approach  not  only  provides  smaller 
expected  values  for  all  MOP' s,  which  means  that  the  results  are  more  accurate,  but  also 
provides  smaller  standard  deviations,  so  the  results  are  more  consistent  as  well. 


D.  TEST  N°3,  TEST  OF  SENSITIVITY  OF  THE  NEW  PROPOSAL 


This  test  tries  to  study  the  dynamical  behavior  of  the  proposed  method.  In  this 
case,  five  clocks  are  simulated.  For  each  run,  an  absolute  covariance  matrix  R  is  first 
generated.  Then,  the  relative  covariance  matrix  S  is  computed  and  finally  R  is  estimated 
by  solving  the  minimization  problem  given  in  (2.43).  The  variance  terms  corresponding 
to  the  elements  of  the  main  diagonal  of  R  are  simulated  according  the  following 
predefined  pattern: 
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r. .  =  15  +  10-cos(^  —  ?); 

11  1024 

_  _  _  ,  2  *  t  TC  x 

r97  =  20  -  8  •  cos( - +  — ) 

1024  2 


r33  =  30 - 


20 -r 
1024 


r44  =  5  +  25  • 


1024 


r55=  15  +  10-cos(-^^; - — ); 


1024 


t  =  {0, 1,  2,...,  1024} 


(3.4) 


Then,  for  i  ^  j ,  the  rest  of  the  elements  of  R  are  computed  using: 

rij  =  Pj  '  (3-5) 

Where  the  correlation  coefficients  p..  are  fixed  at  0.1. 


The  output  of  the  algorithm  is  the  estimation  of  the  absolute  covariance  matrix 
R .  The  idea  behind  the  simulation  is  to  test  the  sensitivity  of  the  algorithm  in  order  to 
check  if  small  changes  in  the  input  parameters  might  produce  big  changes  in  the  output. 
The  total  number  of  runs  in  this  simulation  is  fixed  to  210 


The  usual  approach  to  select  the  weight  for  each  clock  is  to  make  the  weights 
inversely  proportional  to  the  computed  variance  of  each  clock,  as  shown  in  section  II  B. 
In  Figures  7  and  8,  the  simulated  and  computed  variance  terms  (diagonal  elements  of  the 
covariance  matrix)  and  the  corresponding  weights  are  shown.  The  true  value  used  in  the 
simulation  has  also  been  drawn  in  a  solid  line. 
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Simulated  vs.  Computed  Allan  Variances 


run  number 

Figure  7.  Simulated  vs.  Computed  Allan  variances 


Simulated  vs.  Computed  Weights 


Figure  8.  Simulated  vs.  Computed  Optimal  Weights 
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In  all  cases,  the  algorithm  converged  to  a  solution  close  to  the  true  value  and  no 
extra  runs  with  different  initial  conditions  are  needed.  The  mean  value  and  its  standard 
deviation  for  the  relative  error  of  both  Allan  variances  and  weights  are  presented  in 
Tables  14  and  15 


Allan  Var 

Clock  #1 

Clock  #2 

Clock  #3 

Clock  #4 

Clock  #5 

Mean  rel.  err. 

+0.0287 

+0.0255 

+0.0256 

+0.0310 

+0.0304 

Std  rel.  err. 

+0.0434 

+0.0318 

+0.0327 

+0.0478 

+  0 . 0424 

Table  14.  Summary  of  relative  errors  in  Allan  variance 


Weights 

Clock  #1 

Clock  #2 

Clock  #3 

Clock  #4 

Clock  #5 

Mean  rel.  err. 

+0.0022 

+0.0060 

+0.0059 

-0.0005 

+0.0003 

Std  rel.  err. 

+0.0419 

+0.0390 

+0.0386 

+0.0453 

+0.0461 

Table  15.  ! 

Summary  of  rel 

alive  errors  in  Weights 

The  histograms  showing  the  distribution  of  relative  errors,  grouped  by  clock,  have 
been  included  in  Figure  9. 
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Rel  ati  ve  error  range 


Rel  ati  ve  error  range 


Figure  9.  Relative  errors  in  Allan  variances  and  Weights 


The  results  are  clearly  acceptable.  In  77%  of  the  runs  the  errors  are  smaller  than 
0.05%,  and  in  98%  of  the  runs  the  errors  are  smaller  than  0.10% 


One  final  note  about  the  algorithm  is  that,  in  general,  it  seems  to  underdetermine 
the  Allan  variances  associated  with  each  clock.  However,  since  this  effect  happens  with 
all  of  the  clocks,  the  final  result  is  that  for  the  calculations  of  the  weights  the  effect  is 
compensated.  It  can  be  seen  that  in  Figure  7,  the  computed  variances  are  below  the  true 
values  identified  by  the  solid  line.  However,  in  Figure  8  the  computed  weights  are  more 
centered  near  the  true  values,  which  is  beneficial  since  they  are  the  values  needed  to 
characterize  each  clock. 
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E. 


TEST  N°4,  PHASE  SIMULATION 


In  this  case,  a  complete  simulation  is  carried  out  to  analyze  the  whole  process. 
The  differences  in  time  among  the  different  pairs  of  clocks  are  generated  in  advance  in 
order  to  be  later  used  as  the  input  data  for  the  program.  The  process  of  generating  the 
simulated  data  involves  some  decisions  about  the  kind  of  noise  that  is  going  to  be 
introduced  in  the  system  and  which  requires  a  close  study.  Since  this  simulation  involves 
several  stages,  all  of  them  need  to  be  checked.  In  particular,  whenever  intermediate 
results  are  obtained,  they  are  compared  to  the  true  values  of  the  same  parameters.  The 
steps  that  will  be  analyzed  and  checked  in  detail  are: 

•  Generation  of  the  appropriate  type  of  noise  process. 

•  Generation  of  the  fractional  frequency  deviations  and  time  differences 
between  pairs  of  clocks. 

•  Estimation  of  the  relative  covariance  matrix  S  associated  with  the 
ensemble. 

•  Estimation  of  the  absolute  covariance  matrix  R  associated  with  the 
ensemble. 

•  Estimation  of  the  weight  associated  with  each  clock. 

1.  Generation  of  the  Appropriate  Type  of  Noise  Process. 


Atomic  clocks  may  show,  in  general,  the  five  different  types  of  noise  presented  in 
chapter  I.  They  are  characterized  by  the  spectral  components  present,  that  is  by  the 
Power  Spectral  Density  (PSD).  The  different  types  are  easily  identified  in  the  frequency 
domain  by  the  slope  of  the  PSD  of  the  fractional  frequency  deviation  Sv(f)  in  a 

logarithmic  plot.  Fortunately,  there  is  also  a  correspondence  between  the  frequency 

domain  and  the  time  domain.  This  correspondence  implies  that  the  type  of  noise  can  also 
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be  identified  in  the  time  domain  by  plotting  the  Allan  variance  cr(r)  (or  Allan  deviation, 
crv(r))  or  Modified  Allan  variance  Mod  cr(r))  as  a  function  of  the  averaging  time  r  . 


The  first  step  in  this  simulation  is  related  to  the  generation  of  time  data  corrupted 
with  additive  noise  of  specific  spectral  characteristics.  Table  16  summarizes  the  5  types 
of  noise  and  its  characterization  in  both  domains,  time  and  frequency,  and  includes  the 
results  shown  in  Tables  1  and  8. 


Noise  Process 

Slope  Characteristics  in  Log-Log  Plot 

Frequency 

Domain 

Time 

Domain 

sy(f) 

W) 

2 

M  <7y 

Slope 

a 

p 

A 

/ul  2 

A’ 

Random  Walk  Frequency  Modulation 

-2 

-4 

1 

y2 

1/2 

Flicker  Frequency  Modulation 

-1 

-3 

0 

0 

0 

White  Frequency  Modulation 

0 

-2 

-1 

-1/2 

-1/2 

Flicker  Phase  Modulation 

1 

-1 

-2 

-1 

-1 

White  Phase  Modulation 

2 

0 

-2 

-1 

-1/2 

Table  16.  Noise  processes  used 

in  mode 

ling  frequency  instability 

The  five  types  of  noise  processes  are  first  generated  using  the  algorithm  suggested 
in  [29].  The  noise  is  nonnalized  so  that  in  all  cases  the  Allan  variance  corresponding  to 
an  averaging  time  of  r  =  1  is  equal  to  1 . 


In  order  to  test  if  the  generated  noise  agrees  with  the  desired  characteristics,  a 
straight  line  is  fitted  to  the  log-log  plot  of  both  the  power  spectral  density  of  the  fractional 


60 


frequency  deviations  S  (f)  and  to  the  Allan  deviation  <ry(r)  as  a  function  of  the 

averaging  factor.  The  results  for  the  theoretical  and  computed  slopes  and  the 
corresponding  Figures  1 1  and  12  are  also  shown  below. 


Noise  Type 

Freq.  Domain 

Time  Domain 

Theoretical 

Computed 

Theoretical 

Computed 

RWFM 

-2 

-1.80 

1/2 

0.48 

FFM 

-1 

-1.01 

0 

-0.01 

WFM 

0 

-0.00 

-1/2 

-0.51 

FPM 

+  1 

0.78 

-1 

-0.88 

WPM 

+2 

1.79 

-1 

-1.00 

Table  17.  T 

reoretical  and  computed  slopes  for  t 

le  different  noise  processes 

Allan  Deviation  plot  for  different  noise  processes 


Figure  10.  Simulated  noise  processes  and  corresponding  Allan  deviation  plot 


61 


Power  Spectral  Density  of  Fractional  Frequency  S  (f)  for  different  noise  process 


cc 


Log  digital  frequency  (0,1/2) 

Figure  1 1 .  Simulated  noise  processes  and  corresponding  Power  Spectral  Density  plot 


For  Cesium  Time  Standards  located  in  The  Royal  Observatory  of  the  Spanish 
Navy,  the  expected  type  of  noise  for  averaging  time  in  the  order  of  weeks  or  greater  is 
Random  Walk  Frequency  Modulation[6].  Consequently,  this  is  the  type  of  noise  that  is 
used  in  the  simulation  for  all  clocks. 


2,  Generation  of  the  Fractional  Frequency  Deviations  and  Time 
Differences  Between  Pairs  of  Clocks 


Once  the  proper  type  of  noise  has  been  generated  for  each  clock,  the  next  step 
consists  of  the  generation  of  the  absolute  fractional  frequencies  and  time  differences 
between  each  single  time  standard  and  an  ideal  clock.  Once  these  fractional  frequencies 
are  known,  the  absolute  time  differences  between  each  clock  and  an  ideal  clock  can  be 
computed  by  numerical  integration  of  the  fractional  frequencies.  The  last  step  to  finally 
compute  the  relative  time  differences  between  each  pair  of  real  clocks  is  done  by  pairwise 
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subtraction  of  the  absolute  time  differences  previously  obtained.  These  relative  time 
differences  among  the  real  clocks  will  be  the  required  input  to  the  program. 


The  absolute  and  relative  fractional  frequency  fluctuations  and  time  differences 
associated  with  the  ensemble  are  shown  in  Figures  12  to  15.  Since  the  time  origin  does 
not  affect  stability,  all  of  the  clocks  are  initially  supposed  to  be  perfectly  synchronized. 

Absolute  Fractional  Frequency  fluctuations  Clock  #/  -  Ideal  Clock,  /'=1,2,3,4,f 


Figure  12.  Absolute  Fractional  Frequency  fluctuations 
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Absolute  Time  difference  Clock#/'-  Ideal  Clock,  /= 1,2, 3, 4, 5 


Figure  13.  Absolute  Time  differences 


Relative  Fractional  frequency  fluctuation  Clock  #/  -  Clock  #5,  /=1 ,2,3,4 


Sampl  e  poi  nts 


Figure  14.  Relative  Fractional  Frequency  fluctuations 
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Relative  Time  difference  Clock  #/  -  Clock  #5,  /=1 ,2,3,4 


Sample  points 


Figure  15.  Relative  Time  differences 


The  data  shown  in  Figure  12  is  the  input  data  required  to  characterize  each  clock; 
it  represents  how  the  time  scales  provided  by  each  clock  diverge  with  time.  The 
laboratory  data  available  for  this  study  is  similar  to  the  data  shown,  but  coming  from  real 
atomic  clocks  instead  of  simulated. 


3,  Estimation  of  the  Relative  Covariance  Matrix  S  Associated  with  the 
Ensemble 


Once  the  time  series  containing  the  time  differences  between  each  of  the  clocks 
and  the  reference  is  known,  the  relative  Allan  (co)variance  matrix  S  needs  to  be 
computed.  In  order  to  do  so,  some  decisions  need  to  be  made  in  advance  because  there  is 
a  trade  off  between  the  number  of  sample  points  to  use  in  each  computation  of  the  matrix 
S  and  the  desired  accuracy  in  the  time  location  of  the  changes  in  stability  associated  with 
each  clock.  If  the  number  of  sample  points  used  in  the  computation  of  S  is  too  big,  then 
the  matrix  S  will  present  a  smooth  variation  with  time  and  changes  in  the  stability  of  a 
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particular  clock  may  be  either  filtered  out  or  may  not  be  detected.  On  the  other  hand,  if 
the  number  of  sample  points  used  in  the  computation  of  S  is  too  small,  the  matrix  S  will 
be  very  noisy,  and  consequently  the  later  computation  of  the  absolute  covariance  matrix 
R  will  be  affected. 


In  this  simulation,  several  choices  in  the  sample  size  are  used  to  test  this  effect,  in 
particular,  sample  sizes  of  100,  200  and  500  data  points  are  used  to  compute  S  .  The 
effect  on  sample  size  is  shown  in  Figure  16,  along  with  the  theoretical  or  true  value  for 
the  diagonal  and  off-diagonal  terms  of  the  matrix  S.  As  expected,  the  oscillations  of  the 
computed  value  around  the  true  one  decrease  in  amplitude  with  the  number  of  data  points 
used  in  the  computation. 


Simulated  vs.  computed  diagonal  terms  of  matrix  S 


m,  dti.  i  7T  |  jn  i  i  i  i 

Aif-t  1  .  fll  ||  ill  |  |  j  | 

-  sn  : 
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—  833  t 
-  S44 

i,r'  \  l .V  ?f : 

Sample  Sze=  100  1  1 

°0  1000  2000  3000  4000  5000  6000  7000  8000 


Figure  16.  Diagonal  terms  of  matrix  S 

In  order  to  solve  the  minimization  problem  to  find  the  diagonal  terms  of  the 
matrix  R ,  the  complete  S  matrix  must  be  first  computed.  The  off-diagonal  terms  of  S  are 
plotted  in  Figure  17. 
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Simulated  vs.  computed  off-diagonal  terms  of  matrix  S 


Data  Samples 


Figure  17.  Off-diagonal  terms  of  matrix  S 


It  is  no  surprise  that  all  of  the  off-diagonal  terms  of  the  matrix  S  are  very  similar. 
The  reason  is  that,  for  the  case  of  five  clocks,  all  of  them  are  formed  as  a  linear 
combination  of  off-diagonal  terms  of  the  matrix  R  (which  are  small)  plus  the  variance 
tenn  corresponding  to  the  clock  that  is  chosen  as  reference  (r55),  which  can  be  several 
orders  of  magnitude  greater  than  the  off-diagonal  terms  of  R. 


The  results  for  the  matrix  S  shown  in  Figures  13  and  14  are  the  input  values  for 
the  algorithm  that  obtains  the  diagonal  terms  of  R.  Noise  enters  into  the  estimation 
process  of  S,  and  is  consequently  carried  to  the  end  of  the  process.  This  is  the  reason  why 
accurate  values  for  R  are  difficult  to  achieve  in  practice. 
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4. 


Estimation  of  the  Absolute  Covariance  matrix  R  Associated  with  the 
Ensemble 


The  last  step  in  the  simulation  corresponds  to  the  minimization  problem  that 
allows  the  computation  of  the  variance  terms  associated  with  each  clock.  This  process  is 
carried  out  in  two  steps  to  help  the  convergence  of  the  method.  In  the  first  step  the  three 
comered-hat  method  is  used  to  find  an  approximate  solution.  This  solution  is  passed  to 
the  method  to  be  used  as  an  initial  solution  in  the  algorithm  described  in  chapter  II. 


The  final  results,  again,  depend  on  the  sample  size  used.  In  general,  there  is  a 
good  agreement  between  the  simulated  values  for  the  variance  terms  and  the  computed 
values  as  is  shown  in  Figure  18. 


Simulated  vs.  computed  diagonal  terms  of  matrix  R 


Figure  18.  Diagonal  terms  of  matrix  R 
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5. 


Estimation  of  the  Weights  Associated  with  Each  Clock. 


This  is  a  straightforward  step,  since  the  optimal  value  for  the  weights  associated 
with  each  clock  is  proportional  to  the  inverse  of  its  variance.,  After  some  manipulation  of 
the  data  plotted  in  Figure  15,  the  weights  are  computed  and  plotted  in  Figure  16. 


3  mu  I  at  ed  vs.  computed  weights 


Figure  19.  Simulated  vs.  computed  weights 


Figure  19  shows  that  there  is  a  trade  off  between  the  number  of  samples  used  in 
the  computation  and  the  time  resolution  obtained.  If  the  number  of  samples  is  too  small, 
then  the  result  is  noisy.  On  the  other  hand,  if  the  number  of  samples  is  too  big,  time 
changes  in  the  behavior  of  the  clocks  are  not  immediately  detected. 
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IV.  APPLICATION  OF  THE  METHOD  TO  REAL  DATA 


A.  OVERVIEW 


In  this  chapter,  a  data  file  obtained  from  the  Royal  Observatory  of  the  Spanish 
Navy  is  analyzed.  The  file  consists  of  differences  in  time  between  4  atomic  clocks  and 
another  atomic  clock  of  the  same  characteristics  chosen  as  a  reference.  The  measurement 
noise  is  negligible  compared  to  the  noise  introduced  by  the  atomic  clocks,  so  only  4 
columns  of  data  are  required  to  perform  the  study.  For  the  rest  of  the  chapter,  the  clocks 
will  be  referred  to  as  Clock  #1,  Clock  #2,  ...  and  Clock  #5.  The  data  file  has  the 
following  structure,  where  the  differences  in  time  are  expressed  in  nanoseconds: 


Day 

Clock  #1- Clock  #5 

Clock  #2-  Clock  #5 

Clock  #3-  Clock  #5 

Clock  #4-  Clock  #5 

0 

5 

1245 

Tal 

ale  18.  Structure  ol 

~  the  data  file 

B.  INITIAL  EXPLORATION  OF  THE  DATA 


The  raw  data  contained  in  the  file  is  shown  in  Figure  20.  The  parabolic  shape  of 
the  data  show  the  typical  effect  present  in  Cesium  Standards  of  a  constant  linear  drift  in 
frequency,  which  is  equivalent  to  a  parabolic  drift  in  time.  These  are  systematic  effects 
that  must  be  estimated  and  subtracted  from  the  data. 
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Figure  20.  Raw  Time  differences  between  clocks 


C.  INITIAL  FILTERING 


Before  using  the  data,  it  is  necessary  to  preprocess  it  in  two  ways.  First,  the  initial 
time  difference  between  clocks  is  just  a  constant  term  that  plays  no  role  in  the  stability 
characteristics  of  each  clock,  so  it  is  subtracted  so  that  all  of  the  clocks  can  be  considered 
to  be  perfectly  synchronized  in  time  at  the  beginning  of  the  analysis.  Second,  the 
deterministic  effects  associated  with  a  constant  drift  in  frequency  can  also  be  filtered  by 
fitting  a  second  order  polynomial  to  the  data  and  subtracting  it  from  the  raw  data.  When 
both  processes  are  performed,  the  noise  processes  appear  evident.  The  time  residuals  and 
associated  fractional  frequency  deviations  are  shown  in  Figures  21  and  22;  they  are  the 
data  necessary  to  carry  out  the  subsequent  analysis. 


72 


Fractional  frequency  (x  1 0'y/(5  x  86400)  Time  difference  (nanoseconds) 


Filtered  Time  differences:  Clock  #/  -  Clock  #5,  /={  1 ,2,3,4} 


Figure  2 1 .  Filtered  Time  differences  between  clocks 

Filtered  Fractional  frequency:  Clock  #/'  -  Clock  #5,  /={  1 ,2,3,4} 
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Figure  22.  Filtered  Fractional  frequency  between  clocks 
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D. 


CHARACTERIZATION  OF  THE  ENSEMBLE 


The  final  step  consists  of  the  application  of  the  new  method  that  has  been 
developed  in  this  thesis.  The  input  data  are  the  data  shown  in  Figure  19.  In  order  to 
capture  non  stationary  changes  in  the  stability  characteristics  of  the  clocks,  a  window 
length  must  be  chosen  to  compute  the  Allan  variances  of  each  clock  for  the  period  of  time 
corresponding  to  the  window.  For  the  data  provided,  a  window  length  of  128  samples  is  a 
good  trade  off  between  time  accuracy  and  convergence  of  the  method. 


Figure  23  provides  the  computed  Allan  deviation  for  the  data  file  provided  for  this 
work  for  an  averaging  time  of  r  =  5  days. 


Computed  Allan  Deviation,  x=5  days 


Figure  23.  Computed  Allan  deviation  for  the  ensemble 
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For  comparison  purposes,  the  same  available  data  file  of  time  differences  among 
the  clocks  of  the  ensemble  is  also  analyzed  using  a  method  based  on  the  classical  Three- 
Cornered  Hat  approach.  For  this  particular  case,  no  negative  variances  are  found  and 
consequently,  the  results  are  close  to  the  ones  obtained  with  the  approach  presented  in 
this  thesis  However,  when  the  window  width  is  changed  to  64  samples  instead  of  128, 
the  Three-Cornered  Hat  approach  provides  no  results  whatsoever  since  some  of  the 
computed  variances  are  negative.  The  results  are  shown  in  Figure  24. 


Computed  Allan  Deviation  (Three  Cornered  Hat),  x=5  days 


Figure  24.  Computed  Allan  deviation  (Three  Cornered  Hat) 
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CONCLUSIONS 


In  this  thesis,  the  theory  surrounding  the  characterization  of  precision  clocks  and 
oscillators  is  introduced  in  order  to  derive  an  algorithm  to  compute  a  measure  of 
goodness  for  the  stability  of  the  time  scale  of  each  clock  in  an  ensemble. 


The  problem  of  assigning  weights  to  each  clock  is  faced  from  an  operations 
research  viewpoint  by  fonnulating  a  nonlinear  optimization  problem  subjected  to 
nonlinear  constraints  that  are  based  on  physical  considerations. 


One  new  algorithm  to  find  optimal  weights  for  an  ensemble  of  atomic  clocks  is 
introduced,  and  three  different  measures  of  performance  are  also  defined  in  order  to  be 
able  to  perfonn  a  comparison  among  the  proposed  new  method  and  four  other  methods 
actually  used  in  time  laboratories. 


A  comparison  shows  that  the  proposed  new  method  provides  better  values  for  the 
MOP’s,  and  also  that  the  variability  in  the  three  MOP' s  is  the  smallest  using  the  new 
method. 
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